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PREFACE 


The first NASA Computational Fluid Dynamics (CFD) Conference was held at 
NASA Ames Research Center on March 7 through 9, 1989. Conference objectives were 
to disseminate CFD results to industry and university researchers, to promote 
synergy among NASA CFD researchers, and to permit feedback from researchers 
outside of NASA concerning issues facing the discipline of CFD. The focus of the 
1989 conference was on the application of CFD technology. As a result of three 
panel sessions held at the 1989 conference, the following conclusions were 
reached : 


• NASA's program was too heavily focused on applications. 

CFD developers need to better understand the needs of the users. 

Industry needs more reliable and cost effective CFD tools. 

Three critical areas of research were identified: 

-- More algorithm research, particularly for Navier-Stokes solvers with 
unstructured grids. 

-- More research on geometric modeling with rapid, accurate, and 
effective surface definition methods. 

-- More research on grid generation stressing faster, more efficient, 
and improved quality with emphasis on reduced set-up time and 
complexity. 

Additional concerns resulting from the panels may be found in NASA Conference 
Publication 10038, Vol. 1, p. xiv. 

The objectives of the second NASA CFD Conference, held at the Ames Research 
Center on March 12 through 14, 1991, were the same as the 1989 conference: to 
disseminate CFD research results to industry and university CFD researchers, to 
promote synergy among NASA CFD researchers, and to permit feedback from 
researchers outside NASA on issues pacing the discipline of CFD. However, the 
focus of the 1991 conference was on the elements that comprise the discipline of 
CFD, rather than on CFD applications. These elements are 

Algorithm development and analysis 

• Algorithms for real gas /combus t ion modeling 

• Algorithms for advanced computer architectures 

Turbulence modeling/flow physics 

Scientific visualisation 

• Grid generation 

Other CFD methodologies 

The intent of the 1991 conference was to emphasise CFD development underway 
at NASA Centers, rather than the applications activity. 
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This publication is a collection of 20 presentations made by members of the 
NASA Lewis Research Center staff; the members of the Institute for Computational 
Mechanics in Propulsion (ICOMP) ; Sverdrup Corporation, Lewis Research Group; and 
visiting academics at Lewis. Presentation abstracts, (published in bound form 
at the conference) as well as viewgraph material used during the conference are 
included in this publication. The intent is to display research underway on the 
fundamentals of CFD at Lewis. 

The conference was sponsored by the Aerodynamics Division, Office of 
Aeronautics, Exploration and Technology (OAET) , NASA Headquarters, Washington, 
D.C. 


This preface is based on an adaptation of that appearing in NASA Conference 
Publication 10038, proceedings of the CFD Conference held at Ames Research 
Center, 1989. 


Louis A. Povinelli 

NASA Lewis Research Center 

April 8, 1991 
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Development of New Flux Splitting Schemes 

Meng-Sing Liou and Christopher J. Steffen, Jr. 

Internal Fluid Mechanics Division 
NASA Lewis Research Center 

Maximizing both accuracy and efficiency has been the primary objective in 
designing a numerical algorithm for computational fluid dynamics (CFD). This is 
especially important for solution of complex 3D systems of Navier-Stokes equations 
which often include turbulence modeling and chemistry effects. Recently, upwind 
schemes have been well received for both their capability of resolving discontinu- 
ities and their sound theoretical basis in characteristic theory for hyperbolic sys- 
tems. With this in mind, we present two new flux splitting techniques for upwind 
differencing. 

The first method is based upon High-Order Polynomials Expansions (HOPE) of 
the mass flux vector. The present splitting results in positive and negative mass-flux 
components that vanish at M — 0. Thus the error in the Van Leer scheme which 
results in the diffusion of the boundary layer is eliminated. We also introduce several 
choices for splitting the pressure and examine their effects on the solution. 

The second new flux splitting is based on the Advection Upwind Splitting 
Method (AUSM for short). In Navier-Stokes calculations, the diffusion error present 
in Van Leer’s flux-splitting scheme corrupts the velocity vector near the wall. In the 
AUSM, a proper splitting of the advective velocity component leads to an accurate 
resolution of the interface fluxes. The interface velocity is defined using the Mach 
number polynomial expansion in the mass flux, then the convective fluxes follow di- 
rectly. Again, several choices of pressure splitting are possible among which a simple 
Mach number splitting according to characteristics appears to be the best in terms 
of accuracy. The scheme has yielded results whose accuracy rivals, and in some 
cases surpasses that of Roe’s method, at reduced complexity and computational 
effort. 

The calculation of the hypersonic conical flow demonstrates the accuracy of 
the splittings in resolving the flow in the presence of strong gradients. The second 
series of tests involving the 2D inviscid flow over a NACA 0012 airfoil demonstrate 
the ability of the AUSM to resolve the shock discontinuity at transonic speed and 
the level of entropy generation at the stagnation point. 

In the third case we calculate a series of supersonic flows over a circular cylin- 
der. The Roe splitting in all conditions and grids tested yielded anomalous solu- 
tions (sometimes referred to as the carbuncle phenomenon), which could appear as 
nonsymmetric, protuberant, or indented contours, see Figs. 1. The AUSM gave 
expected solutions in all calculations. 

The fourth test deals with a 2D shock wave/boundary layer interaction. This 
provides an opportunity to accurately resolve a laminar separation region and to 
compare the ability to resolve a non grid-aligned shock with other methods, as 
shown in Figs. 2. 
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Fig. 1 Mach contours for a Mach 6 flow over a circular cylinder: (a) the AUSM 
solution, and (b) the Roe solution displaying a protuberant, two-shock solution. 
Note the calculations were performed for the front half plane, thus no symmetry 
was imposed. 
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Fig. 2 Mach contours for a shock wave/boundary layer interaction problem: (a) 
the AUSM solution showing sharp resolutions of various waves, and (b) the Roe 
solution showing more diffusive waves. 
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Motivation 

To couple Accuracy and Efficiency together within an Upwind 
Flux-Splitting Scheme 


Two Possible Improvements 


1. Less Diffusive FVS Technique 

• High Order Polynomial Expansion (HOPE) 

• Unresolved Stability Concerns 

2. New Adveclive Velocity Splitting Technique 

• Advection Upwind Splitting Method (AUSM) 

• Simple, Accurate, Efficient and Robust 
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Quasi-2D Viscous Conical Flow 




Adiabatic Wall Boundary Condition 
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Euler Equations 
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Conical Flow: First Order Accurate Results 


Pressure Profile 



Temperature Profile 


15 


M 


1:1 


12 


1 1 


10 


Van 1 .eer 
Roe 

. AUSM 
M - 7.95 
Pr =■= 1.0 


1 %. 


a 

T/T. 


12 


16 


») theory 


23 





Conical Flow: First Order Accurate Results 


Convective Velocity Profile 



* Van Leer 
o Roe 

• AUSM 

M = 7.95 
Pr = 1.0 


1.030 - 0.020 - 0.010 0.000 

u* A. 


o.oio 


2D Shock Wave - 
Boundary Layer Interaction 

Geometry 

Impinging Shock Reflecting Shock 



Re = 2.96 x 10 5 
Shock Angle = 32.58° 
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norm 


Conical Flow: Convergence Histories 


First Order Accurate Scheme Second Order Accurate Scheme 



iterations iterations 

Shock Wave/Boundary Layer: Mach Number Contours 


AUSM Scheme 0(A 2 ) Roe FDS Scheme 0(A 2 ) 
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Mach 6.0 Blunt Body Flow: Carbuncle Phenomenon 
Velocity Vectors 



Mach 3.0 Blunt Body Flow: Mach Number Contours 


AUSM Scheme 0(A) Grid Roe FDS Scheme 0(A) 
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Supersonic Blunt Body Flow 



Mach 6.0 Blunt Body Flow: Mach Number Contours 
Roe FDS Scheme 0(A) Grid AUSM Scheme 0(A) 
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Shock Wave/Boundary Layer: 


Friction Coefficient 



Pressure at the Wall 
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♦Experiment by Hakkinen, Greber, Trilling and Abarbanel 


Conclusions 


• Accuracy of AUSM rivals FDS schemes (i.e. Roe Splitting) 


• Efficiency and Simplicity of AUSM rivals FVS schemes 
(i.e. VanLeer Splitting) 

• AUSM is a new family of Upwind Flux Splitting Tech- 
niques 
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AN ITERATIVE IMPLICIT DDADI ALGORITHM 
FOR SOLVING THE NAVIER-STOKES EQUATION 


S.C. CHEN, N.S. LIU, and H.D. KIM 
INTERNAL FLUID MECHANICS DIVISION 
NASA LEWIS RESEARCH CENTER 


An algorithm utilizing a first-order upwind split flux 
technique and the diagonally dominant treatment is proposed to 
be the temporal operator for solving the Navier-Stokes 
equations. During the factorization process, if care is taken 
to ensure the symmetry of the operator, the resulting 
algorithm will be stable when spatial derivatives of the right 
hand side residual are evaluated by the central differencing 
scheme as well as the upwind differencing Roe scheme. 

Temporal accuracy is assured by the inner Newton-iteration 
process first introduced by Rai (1986) . 

Given the limit of a five-point stencil, the right hand side 
flux derivatives are formulated by several commonly used 
central and upwind schemes. Their performances are studied 
through a test case of free vortex convection in a uniform 
stream. From these results, a superior treatment for 
evaluating the flux term is proposed and compared with the 
rest. 

The application of the proposed algorithm to the full Navier- 
Stokes equations is demonstrated through a calculation of flow 
over a backward facing step. Results are compared against the 
calculation done by using the fourth-order central 
differencing scheme with artificial damping. 


29 



LEWIS RESEARCH CENTER 


OBJECTIVE 

GIVEN A STENCIL OF 5 POINTS, FIND A 
SCHEME THAT COMBINES THE MERITS 
OF BOTH CENTRAL DIFFERENCING AND 
UPWIND DIFFERENCING SCHEMES. 
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Algorithm: In 2 Consecutive Steps 

* (n) 

[a 1+ DAt + S At + (LAt)j AQ = At • RHS 


(n) 

[a I+DAt + SAt] AQ + [At L(AQ) + At U (AQ)] 

(n) 

= At - RHS 

*. (n+1) 

[a 1 + D At + S At + (u At)] AQ = At-RHS 


(n + 1) ** ** 

[al+ D At+S At] AQ + [At L (AQ) + At U(AQ)] 

(n+1) 

= At - RHS 
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Formal Error in 2 Time Steps is 


2'Err— A t'-U 


^5AQ A 
V 5t J 


+At 2 ■ {L [a I + At(D + S)] _1 U} (AQ) 


-Ar-L 


'saq' 

v 5t y 


-1 

At 2 • *{ll [a I + At (D + S )] Lj* (AQ) 


Define 



(n) 

AQ 

At 


(a-1) 


(n-1) 

AQ 

At 


P 

= a X 

. q=i 


(q) 

AQ 

At 


(a-1) 


(n-1) 

AQ 

At 


Then 


RHS = 


5Q 

5t 


(P) 


a 


AQ 

At 


6E_5F_5G 5(Ev) 8(Fv) 8(Gv) 
5^ 8n K + 5^ + Sri + 8 £ 
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Geometry 


■< 


Ax = Ay = 1/4 
At = 0.05 


(Vortex core radius = 1) 

f P7 = 0.72 


CFL ~ 0.7 
M°© = 0.5 

r = -0.5 

Innerm = 4; NORDR = 2 
Invisid 


7 = 1.4 
Too = 519 °R 

p = 2.374 x 1 0" 3 slug/ft 3 

OO 

Air 



U,V, P, T = f„(t) 

Travel a total of 45 length in 900 steps. 


Test 1 4th Order Central 

Flux derivative in 1-D difference for 



2_ 
3 h 


( E i+r E i-i) 


12 h 


( E i+2 E i - 2 ) 


N ,2 j (4 — 2 n ~ 1 > 

~ r 3 n! 
rr=5 0 


Ef n) h n-1 (N-> 00) 

1 ' 


* Contains only the odd derivatives leading error is of 
order "4" 
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Pressure 






Consider: 


E'j = [4th order centra! differencing] 


N.2 1 (4 - 2 n_1 ) 


"X q 

n = 5 3 


E (") h "- 1 
ni 1 


and 


E'j = [3rd order upwind biased differencing] 


(A) 


N.2 1 (4-2 n_1 ) 


1 T 
n = 5 3 


E (") h n-1 
ni i 


N.2 -i (2-2 11 - 1 ) („) hn _! 


I o 
n = 4 3 


m 


(B) 


The Dissipative Equivalence is 

Diss = (A) - (B) 

Diss = - N i 2 i (2 ~ n 2 — * E^ n) h"- 1 

= [3rd order upwind difference] 

- [4th order central difference] 
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8E 

Let — = 4 th central + y • DISS 
5x ' 

Then = 

Ej = {4th order central differencing 

r N,2 1 (o~o n ~^) 

+y ■ l- I l — - E<") h"- 1 ]) 
n= 4 3 n! i 

_N ,2 1 ( 4 _ 2 n ~1 ) (n) n _t 

n = 5 3 n! 

= {(1- y)‘ [4th order central difference] 

+ y • [3rd order upwind biased difference]} 

+ Trunc. error 

Test 5. 3rd Order Upwind Biased Scheme 

E i = [ih (E w- E w ) + ii; <- 3E i +4E i + r E 42). 

N ^ 2 1 (4-2 n_1 ) i n) .n-1 
- X “ — — E. h 

n= 5 3 ■ 1 (odd derivatives) 

_N,Z 1 (2-2 n -1) 

^ 3 ni i 

n=4 (even derivatives) 

* Leading error is of order "3" 

Test 6 

Set 7 = 0.1 8 so that leading error term has a coefficient of 0.01 5 


38 



Vorticity Magnitude 
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Case 2 


Flow over a backward-facing step 


Re d 100 100 

389 389 

500 1000 

600 
1000 

t t 

4th order central Proposed 

+ artificial viscosity scheme 

0 = 0.5) (7*0.5) 


Max. L 2 - Residual < 5x1 O ' 5 for all cases 

f Moo =0.1 8 CFL = 50.0 

J Pr = 0.72 NORDR = 1 

j 7 = 1 .4 Innerm = 1 

^ Air; viscous flow 


U, T given 
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Reattachment and Detachment Length I 



Ref.-JFM: voll 27, pp. 473-496. ArmaJyetaL 


Reattachment and Detachment Length II 


25 r 


20 


Solid = Proposed scheme 
Hollow = 4th central 


R 2 


5 15 

o 

LJ 

X 

£L 

U1 

in 

II 

O 10 


o 

2 


& 


9 R i 

« S 2 


4 


l I t 1 1 L 

0 200 400 600 800 1000 t200 
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□ Measurements 

Proposed scheme 

- - - 4th Order central 


V - Velocity profiles 
Re = 100 



□ Measurements 

Proposed scheme 

4th Order central 


V - Velocity profiles 
Re = 389 
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Static Pressure - Proposed Scheme 



V Component - 4th Order Central 




-3 


0 


Re = 389 
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a 


37 


Re = 600 
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V Component - Proposed Scheme 




(7= NASA 


"" """ LEWIS RESEARCH CENTER 

CONCLUDING REMARKS 


* CENTRAL DIFFERENCING: 


SUSCEPTIBLE TO SPURIOUS OSCILLATIONS, 
ARTIFICIAL DAMPING NEEDS TO BE ADDED. 

* UPWIND DIFFERENCING: 


INTRINSIC YET EXCESSIVE NUMERICAL DISSIPATION. 

* PROPOSED SCHEME: 

USES 4th ORDER CENTRAL DIFFERENCING SCHEME WITH 
INTRINSIC NUMERICAL DISSIPATION PROVIDED BY THE 
UPWIND DIFFERENCING SCHEME. 


INTERNAL FLUID MECHANICS DIVISION 
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The Proteus Navier-Stokes Code 

Charles E. Towne and John R. Schwab 
Internal Fluid Mechanics Division 
NASA Lewis Research Center 



An effort is currently underway at NASA Lewis to develop two- and three-dimensional Navier-Stokes codes, called 
Proteus, for aerospace propulsion applications. The emphasis in the development of Proteus is not algorithm 
development or research on numerical methods, but rather on the development of the code itself. The objective is to 
develop codes that are user-oriented, easily-modified, and well-documented. Well-proven, state-of-the-art solution 
algorithms are being used. Code readability, documentation (both internal and external), and validation are being 
emphasized. 

Proteus solves the Reynolds-averaged, unsteady, compressible Navier-Stokes equations in strong conservation law 
form. Turbulence is modeled using a Baldwin-Lomax based algebraic eddy viscosity model (AIAA Paper 78-257.) 
The governing equations are written in Cartesian coordinates and transformed into generalized nonorthogonal 
body-fitted coordinates. They are solved by marching in time using a fully-coupled ADI solution procedure (Briley 
and McDonald, J. Comp. Phys., Aug. 1977) with generalized first- or second-order time differencing (Beam and 
Warming, AIAA J., Apr. 1978.) The boundary conditions are also treated implicitly, and may be steady or 
unsteady. All terms, including the diffusion terms, are linearized using second order Taylor series expansions. 

In addition to solving the full set of Reynolds-averaged equations, options are available to solve the thin-layer or 
Euler equations, and to eliminate the energy equation by assuming constant stagnation enthalpy. Artificial viscosity 
is used to minimize the odd-even decoupling resulting from the use of central spatial differencing for the convective 
terms, and to control pre- and post-shock oscillations in supersonic flow. Two artificial viscosity models are avail- 
able - a combination implicit/explicit constant coefficient model (Stegcr, AIAA J., July 1978), and an explicit non- 
linear coefficient model designed specifically for flows with shock waves (Jameson, Schmidt, and Turkel, AIAA 
Paper 81-1259.) At NASA Lewis, the code is run on the Cray X-MP and Y-MP computers, and is highly vector- 
ized. 

An extensive series of validation cases have been run, primarily using the two-dimensional planar/axisymmetric 
version of the code. Several flows were computed for which exact solutions to the Navier-Stokes equations exist, 
including fully-developed channel and pipe flow, Couette flow with and without a pressure gradient, unsteady 
Couette flow formation, flow near a suddenly accelerated flat plate, flow between concentric rotating cylinders, and 
flow near a rotating disk. Additional validation cases that have been successfully run include flat plate laminar and 
turbulent boundary layers, boundary layers with confined separation, 2-D and 3-D driven cavities, normal and 
oblique shock-boundary layer interactions, steady and unsteady flows past a circular cylinder, developing laminar 
and turbulent pipe flows, and steady and unsteady flows in a transonic diffuser. The figure shows computed Mach 
number contours and the static pressure distribution for flow through a transonic diffuser. 

The two-dimensional version of Proteus has recently been released, and the three-dimensional code is scheduled for 
release in late 1991. The documentation consists of a three-volume user’s manual (Towne, Schwab, Benson, and 
Suresh, NASA TM’s 102551-102553). Volume 1 is the Analysis Description, and presents the equations and solu- 
tion procedure. It describes in detail the governing equations, the turbulence model, the linearization of the equa- 
tions and boundary conditions, the time and space differencing formulas, the ADI solution procedure, and the 
artificial viscosity models. Volume 2 is the User’s Guide, and contains information needed to run the program. It 
describes the program’s general features, the input and output, the procedure for setting up initial conditions, the 
computer resource requirements, the diagnostic messages that may be generated, the job control language used to 
run the program, and several test cases. Volume 3 is the Programmer’s Reference, and contains detailed informa- 
tion useful when modifying the program. It describes the program structure, the Fortran variables stored in common 
blocks, and the details of each subprogram. 
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LEWIS RESEARCH CENTER 


THE PROTEUS NAVIER-STOKES CODE 


C. TOWNE 

Inlet, Duct & Nozzle Flow Physics Branch 


J. SCHWAB 

Turbomachinery Flow Physics Branch 
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NAVIER-STOKES ANALYSIS 
FOR TRANSONIC DIFFUSER 


Qiwiuii 



"% § 


BOTTOM WALL - PROTEUS 

TOP WALL - PROTEUS 

O BOHOM WALL - EXPERIMENT 

□ TOP WALL - EXPERIMENT 


X— COORDINATE 



PROTEUS NAVIER-STOKES CODE 


• Objective 

— Develop user-oriented, easily-modified, well-documented, 

2- and 3-dimensional Navier-Stokes codes for aerospace 
propulsion applications. 

• Approach 

— Use well-proven, state-of-the-art algorithms. 

— Use a consistent, modular code structure. 

— Emphasize documentation, code readability, and validation. 


PROTEUS NAVIER-STOKES CODE 

ANALYSIS SUMMARY 

• Reynolds-averaged, unsteady, compressible 
Navier-Stokes equations. 

• Strong conservation-law form. 

• Generalized nonorthogonal body-fitted coordinates. 

• Convection and diffusion terms linearized 
using second-order Taylor series expansion. 

• Baldwin-Lomax algebraic turbulence model. 

• Fully-coupled ADI solution procedure, 

Beam-Warming generalized time differencing. 

• Implicit steady/unsteady boundary conditions. 
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PROTEUS NAVIER-STOKES CODE 

CODE FEATURES 


• 2-D, axisymmetric with swirl, or 3-D flow. 

• Thin-layer and Euler options. 

• Wide variety of boundary conditions. 

• Constant stagnation enthalpy option. 

• Constant-coefficient or adaptive artificial viscosity. 

• First- or second-order time differencing. 

• Variety of time step selection methods. 

• Output files for CONTOUR and PLOT3D. 

• Highly vectorized for Cray computers. 

• Extensively commented. 

• Three-volume documentation set. 

— Analysis Description. 

— User’s Guide. 

— Programmer’s Reference. 


PROTEUS NAVIER-STOKES CODE 

VALIDATION CASES 

• Exact Navier-Stokes solutions. 

— Fully-developed channel and pipe flow. 

— Couette flow w/wo pressure gradient. 

— Couette flow formation. 

— Flow near a suddenly-accelerated flat plate. 

— Flow between concentric rotating cylinders. 

— Flow near a rotating disk. 


50 



PROTEUS NAVIER-STOKES CODE 

VALIDATION CASES 

• Blasius flat plate boundary layer. 

• Flat plate boundary layer with confined separation. 

• Turbulent flat plate boundary layer w/wo heat transfer. 

• Flow over a rearward-facing step. 

• 2-D and 3-D driven cavity. 

• Normal and oblique shock formation. 

• Steady and unsteady flow past a circular cylinder. 

• Flow past an airfoil. 

® Developing laminar and turbulent flow in a 
channel and pipe w/wo swirl. 

• Developing 3-D flow in a square duct. 

• Steady and unsteady turbulent flow in a transonic diffuser. 


rr 

Q) 

C 

v_ 

o 

o 

0 

>* 
■ +— « 

J3 

1 
CO 
co 

'co 

as 

CD 



51 




Circumferential Location ft, Degrees 
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STATIC PRESSURE 


NAVIER-STOKES ANALYSIS 
FOR TRANSONIC DIFFUSER 



PROTEUS NAVIER-STOKES CODE 



HISTORY 

• 4-86 

Start. 

• 9-87 

First master file. 

• 9-88 

Emphasis placed on 2-D 
code development. 

• 9-89 

Version 1 .0 of 2-D code released. 

• 3-90 

Documentation published for 
version 1 .0 of 2-D code. 

• 9-90 

Version 1 .2 of 2-D code released. 

• 3-91 

Preliminary Version 1 .0 of 3-D code 
ready for testing. 
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PROTEUS NAVIER-STOKES CODE 
2D Version 



PROTEUS NAVIER-STOKES CODE 
3D Version 
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PROTEUS NAVIER-STOKES CODE 

CURRENT STATUS 

• Version 1 .2 of 2-D code available. 

• 3-D code and documentation being upgraded. 

• Two-equation k-z turbulence model 
available in 2-D code. 

• Grid patching (embedded boundary) capability 
available for 2-D code. 

• Validation studies underway. 

— Laminar/turbulent compressible 
flat plate boundary layers 
— Laminar/turbulent flow past an airfoil. 

— Turbulent flow in transonic diffusers. 


PROTEUS NAVIER-STOKES CODE 

NEAR-TERM PLANS 

• Complete upgrade of 3-D code and documentation. 

• Investigate additional algebraic and pde 
turbulence models. 

• Create diagonalized version. 

• Continue validation of 2-D code, emphasizing 
turbulent flow and flow with heat transfer. 

• Begin more intensive validation of 3-D code. 

• Release 3-D code in fall of 1991 . 


55 



PROTEUS NAVIER-STOKES CODE 

LONG-TERM PLANS 


• Continue validation of 2-D and 3-D codes. 

• Add zonal gridding capability. 

• Investigate convergence acceleration techniques. 

• Add algebraic Reynolds stress turbulence model. 

• Investigate adaptive mesh schemes. 

• Hold user’s workshop after 3-D code release. 


56 




ACCURATE UPWIND-MONOTONE (NONOSCILLATORY) 
METHODS FOR CONSERVATION LAWS 


Hung T. Huynh 

Internal Fluid Mechanics Division 
NASA Lewis Research Center 


The well-known MUSCL scheme of Van Leer (J. Comp. Phys., 1979) — a 
second-order extension of Godunov’s first-order upwind-difference method— is con- 
structed using a piecewise linear approximation. In order that the solutions are free 
from oscillations, a constraint is imposed on the slopes of the piecewise linear recon- 
struction. The MUSCL scheme is second-order accurate at the smooth part of the 
solution except at extrema where the accuracy degenerates to first-order due to the 
monotonicity constraint. Harten and Osher (SIAM J. Numer. Anal., 1987) resolved 
this loss of accuracy by introducing the concept of nonoscillatory reconstruction. 
Their UN02 scheme involves no constraint, is uniformly second-order accurate, and 
does not create oscillations. The UN02 scheme is, however, very diffusive, e.g., it 
smears contact discontinuities. 

To construct accurate schemes which are free from oscillations, we first intro- 
duce the concept of upwind-monotonicity. Next, in a geometric framework in which 
the median function plays a crucial role, we present Van Leer’s constraint and the 
UN02 scheme. In this framework, the constraint is extended in such a manner that 
extrema are not “clipped”, while upwind- monotonicity is preserved. At monotone 
part of the data, the new constraints reduce to that of Van Leer for the second- 
order case, and in the third order case, that of Colella and Woodward (J. Comp. 
Phys., 1984). Several classes of schemes, which are upwind-monotone and of uni- 
form second or third-order accuracy, are then presented. These schemes also satisfy 
the nonoscillatory property of the UN02 scheme. It is shown that the MUSCL, 
UN02 and all the popular TVD (total variation diminishing) schemes (Harten, J. 
Comp. Phys., 1983; or Sweby, SIAM J. Numer. Anal., 1984) are members of the 
new class. Moreover, each TVD scheme corresponds to a uniform second-order ac- 
curate upwind-monotone scheme. The gain of accuracy is obtained by a few extra 
lines of Fortran programming. These new schemes are also extensions of the au- 
thors SONIC schemes (NASA TM 102010, 19S9; also Yokota & Huynh, NASA TM 
102354, 1990) 

Results for advection with constant speed are shown below. The initial profile 
consists of a sin 2 wave, a square wave, a triangle wave and a semi-ellipse wave. 
Each wave contains 20 grid points on a uniform grid of 200 points. Using periodic 
boundary conditions, the profile is advected one period with CFL number 0.5. The 
results after 400 time steps are shown in Figs. (A) for a member of the new class, 
(B) for the UN02, (C) for the PPM (Colella and Woodward, J. Comp. Phys., 
1984), and (D) for Roe’s “Superbee” schemes. It can be seen that the new scheme 
compares favorably with the state-of-the-art methods. 



Upwind Shock-Capturing Methods 

Why new methods? 

TVD (Total Variation Diminishing) schemes 
+ Capture shocks with high resolution 
— First-order accurate near extrema 
UN02 scheme 

+ Uniformly second-order accurate 
+ Nonoscillatory, i.e., no new strict extrema 
— Somewhat diffusive 

New Class of Methods 

New geometric framework 

Existing schemes (TVD, UNO) 

are presented in the new framework 

New schemes 

* Uniformly second-order accurate 

* Upwind-monotone and nonoscillatory 

* Extensions of TVD and UN02 schemes 


du du 
dr + a dx 



a > 0 


u(x, 0) = uq(x). 


Approach 

Projection or Reconstruction, 

Evolution or Upwinding + Time stepping 
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The task for piecewise linear reconstruction: 
given the data {w., }, define (slope) ■ 



The problem: discontinuity 



Van Leer’s constraint 

The linear reconstruction in cell j lies in 
[min(uj_i, Uj, Uj+i), m&x(uj-i,Uj,Uj+i)] 


i i 






a 


m (s, t) = median (s, t, 0) 

= ^[sgn (s) + sgn (£)] min(|s|, \t\) 
to A 0 t 0 

* 1 — 1 X 1 -» 


t 
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The pictures 



Constraint: (slope)^ lies between 0 and 2 r, 

r — m (s, t). 

Algorithm: (slope)- m [ (slope)^ , 2r ] 

1 1 

0 1L\, 

TVD schemes 
Minmod scheme 

(slope) • = r = m(s,t) 

Average limiter 

(slope)^ = m 

Superbee scheme 

(slope)- = m [maxmod(s,i),2r] 

maxmod (s,t) = ^[sgn (s) + sgn (£)] max(|s|,|t|) 

Zj 
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The Peter-Paul Principle 


✓ 



Resolve shock well (Paul), lst-order at extrema (Peter). 
For TVD schemes, u ™ +1 lies between u ” and Uj_ 1 . 

The data are monotone in [xj,Xj+i] if 

Uj ~ i < Uj < Uj+t < Wj+2, Or > Uj > Uj + 1 > Uj+2- 

A scheme is upwind-monotone if monotonicity of the data 
in [xj-\ } Xj] implies u ™ +1 lies between Uj and Uj- 


UN 02 scheme 



(slope) ■ = r = m(s,t) 

+ Second-order accurate 
+ Nonoscillatory 
— Diffusive 
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New methods 


New constraint: 

(slope) lies in I[ 0, 2 r, f] = J[0, r max ] 

r max = sgn(f) max (2|r|, |f|). 

New algorithm: (slope) <— m [ (slope) , r max ]. 

UNO (slope) j = f = m(s, t) 

SONIC-A (slope). = m(^,r raax ) 

SONIC-B (slope)- = m[maxmod(s,t),r max ] 
Second-order accurate and upwind-monotone 



(E>) FROMM SCHEME 
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Conclusions 

A new framework for shock-capturing methods was 
presented 
. Simple 
. Geometric 

. Uses the median function 

The approach leads to new classes of methods: 

. Compare favorably with existing methods 
. Resolve discontinuities with high resolution 
. Uniformly second-order accurate 
. Upwind-monotone (therefore, nonoscillatory) 

. Generalization of TVD and UN02 schemes 
. Vectorize on Cray 
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Numerical Simulation of Conservation Laws 


Sin Chung-Chang 
Internal Fluid Mechanics Division 
NASA Lewis Research Center 

Wai Ming To 

Sverdrup Technology, Inc. 
Lewis Research Center Group 
NASA Lewis Research Center 


A new numerical framework for solving conservation laws is being 
developed. This new approach differs substantially from the well established 
methods, i.e., finite difference, finite volume, finite element and spectral 
methods, in both concept and methodology. The key features of the current 
explicit scheme include: (a) direct discretization of the integral forms of 
conservation laws, (b) treating space and time on the same footing, (c) flux 
conservation in space and time, and (d) unified treatment of the convection and 
diffusion fluxes. 

The model equation considered in the preliminary study is the standard 
1-D unsteady constant-coefficient convection-diffusion equation. In a stability 
study, it is shown that the principal and spurious amplification factors of the 
current scheme, respectively, are structurally similar to those of the 
leapfrog/DuFort-Frankel scheme. As a result, the current scheme has no 
numerical diffusion in the special case of pure convection and is unconditionally 
stable in the special case of pure diffusion . 

Assuming smooth initial data, it will be shown theoretically and 
numerically that, by using an easily determined optimal time-step, the accuracy 
of the current scheme may reach a level which is several orders of magnitude 
higher than that of the MacCormack scheme, with virtually identical operation 
count. 
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Introduction 


The focus of the current method development is entirely on the 
integral form of conservation laws. Key features of the current 
explicit scheme include: 

• Direct discretization of the integral forms of conservation laws. 

• Treating space and time on the same footing. 

• Flux conservation in space and time. 

• Unified treatment of the convection and diffusion fluxes. 

• Low-operation count and high accuracy. 


Unsteady 1 -D Convection-Diffusion Equation 



or ^ ^ 

f S(V) h-ds=0 ; 

V = Any volume in space-time. 

S (V) = Boundary of V. 

— > — > — ) . — > 

h • ds = Flux of h leaving V through surface element ds. 
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Basic Convervation Domains (BCD's 


For each BCD, it is required that 


Sum of incoming fluxes = sum of outgoing fluxes 


As a result, <1 > is also true for the union of any collection of BCD s. 


Space - Time Mesh and BCD's 




dt J along a j mesh line 


x" = jAx + n • bAt, t n = nAt 










Since total flux leaving D'(j,n) = 0 
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Invariant Property (1) 


Let 

u = u G (x,t; a, ft) 
be a solution of 


B 2 u 


Bu Bu 

_ + 3“ — II 

Bx F dx 2 


= 0 


< 1 > 


< 2 > 


Let x* and t* be any real constants. Then 

u = u Q (2x ir -x, 2t,-t; a, -|x) <3> 


tx y t; 


l 

I ' 
i / 

\> 

/ i 

/ i 
/ » 


is also a solution of <2>. As an example, consider 

u Q (x, t; a, p,) — e _4 ‘ n2 l J ' t sin [2-ir(x-at)] <4> 

The above property may be referred to as the p-t-fjt (parity-time 
reversal-p. reversal) invariance. 


Invariant Property (2) 


The main scheme of the Leapfrog/DuFort-Franke! method 
i.e., 


u n+1 _ u n-1 u n 
> J +a- J+1 


u? 


M 


UV 


2At 


2Ax 


= M-- 


j+1 u j 


u n+1_ u n-1 +u n 


J-1 


(Ax)' 


<5> 


is also p-t-p, invariant. In other words, if 

u " = y 0 ( i* n; a >^) 

satisfy <5>, then so does ( j and n^ are constant integers ) 
u" = u 0 (2j-j, 2n- n; a-p,) 

The present scheme is p-t-p. invariant also. Other classical 
schemes generally are not p-t-|x invariant, p-t-u, invariance 
has an impact on the stability property. 
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Stability Condition 


1. The stability condition for the present scheme is 

1 >T 2 

i.e., stability is not dependent on the viscosity fx. 

2. The present scheme has no numerical diffusion in the 
special case of pure convection, and is unconditionally 
stable in the special case of pure diffusion. 


Inconsistency 


• A smooth solution which satisfies the present or the Leapfrog/ 
Dufort-Frankel scheme at the mesh points will satisfy 


3u/dt + a 3u/5x-|x5 2 u/3x 2 = 0 (H*0) <1> 


as Ax -» 0 and At/Ax -» 0. 


• A characteristic line of <1 > is 
t = constant 






Thus a discrete solution may converge to an analytical 
solution only if (At/ Ax) -> 0, 


• For a typical classical scheme, stability requires that 

jtAt jjXAt/Ax) k e bounded. Thus, (At/ Ax) ->0 as A x 0 
(Ax) 2 Ax 
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Problem Definition (21 


A problem is defined by specifying the values of (i) physical 
parameters a and fx, (ii) the mesh parameters b, At and 
Ax (= 1/K), and (iii) total running time t. The total time-step 
number 

def 

n = IN(t/At) — the integer nearest to the ratio (t/At) 

Error of a problem is measured by: 




Accuracy of present method is optimized when either 
(i) x = 0 or (ii)1-t 2 = V3 8.1 t=^ f (a-b)At/Ax, 5 M 4jxAt/(Ax) 2 


„.def , 

E(b,n,At,K) — log 10 


uP - u_(xP,t n ) 
J aV J 



Courant number, r = foo At 
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E(b, 390,1/(45 V3), 30) 



Special Property 

Let K be the number of spatial mesh points. Then the 
present marching scheme may be expressed as 

u n+ 1 = A u n + B n , n = 0,1,2,. .. 

Here (i) u n and u n+1 , respectively, are the 2K x 1 column 
matrices representing 2K numerical variables at the time 
levels n and n+1, (ii) 6 n is a given 2K x 1 column matrix at 
the time level n, and (iii) A is a give 2K x 2K sparse matrix. 


Unusual property: A -1 is also sparse and is known explicity 
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Conclusions (1) 


• Flux conservation is the focus of method development. 

• Unified treatment of space and time, also of convection 
and diffusion fluxes. 

• Stability is not dependent on the viscosity p.. The 
present scheme has no numerical diffusion if §x = 0 
and is unconditionally stable if a = 0. 

• The sparse matrix A in the present marching scheme 
u n+1 = A u n + 6 n has the unusual property that A" 1 is 
also sparse and is known explicity. 

• The present scheme has the same invariant property of 
the physical equation it models. 


Conclusions (2) 

• The present scheme is most accurate and stable when 
t = 0, i.e., b = a 

• For the present scheme, mesh refinement may be made 
locally without impacting global structure. 

• Since initial data required include both u ° and (<9u/<9x)9, 
initial-value specification of the present method is more 
accurate than that of traditional methods. 

• Assuming smooth, initial data, it is shown that, by using 
an easily determined optimal time-step, the accuracy of 
the current scheme may reach a level which is several 
orders of magnitude higher than the MacCormack 
scheme, with virtually identical operation count. 
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A New Lagrangian Method for Real Gases at Supersonic Speed 


C.Y. Loh and Meng-Sing Liou 
Internal Fluid Mechanics Division 
NASA Lewis Research Center 


Abstract 


With the renewed interest in high-speed flights, the real gas effect is of theoretical as 
well as practical importance. In general, this real gas phenomenon is seen in the regions of 
high gradients and discontinuities. The ability of numerically capturing them accurately, 
viz with minimal numerical dissipation and oscillations, has been the main challenge to 
the computational fluid dynamist as well as numerical analyst over the past four decades. 
In the past decade, upwind splittings or Godunov-type Riemann solutions have received 
tremendous attention and as a result significant progress has been made both in the ideal 
and non-ideal gas. However, almost all of these efforts have been exclusively formulated in 
the framework of Eulerian description of fluid motion. Although the Lagrangian descrip- 
tion has been used in several reports, they are always followed by a remap step from a 
Lagrangian grid back to the original Eulerian one. This not only involves tedious interpo- 
lation procedure but also can corrupt the accuracy which was gained by the choice of the 
Lagrangian approach over the Eulerian one. 

However, since all upwind schemes used by today’s major CFD codes are strictly 
formulated in 1-D framework, and only formally extended to multi-dimensions in space. 
Consequently, the attractive property of crisp representation of discontinuities is lost and 
search for genuine multi-dimensional approach has just been undertaken by several leading 
researchers, but still based on Eulerian description. In this paper, we propose a new 
alternative that is formulated using the Lagrangian description, for the calculation of 
supersonic/hypersonic real-gas inviscid flows. This new formulation avoids grid generation 
step which is automatically obtained as solution procedure marches in the “time-like” 
direction. As a result, no remapping is required and the accuracy is faithfully maintained 
in the Lagrangian level. Moreover, the ability of preserving the sharpness of discontinuities 
in muti-dimensions is vividly demonstrated. In this paper, we will give numerical results for 
a variety of real-gas problems consisting of essential elements in the high-speed flows, such 
as shock waves, expansion waves, slip surfaces and their interactions. Finally, calculations 
for flows in a generic inlet and a nozzle are presented. 



Outline 


. New Lagrangian formulation and the conservation form 
for 2-D steady supersonic flow (for perfect gas) 

. Advantages of the new Lagrangian method 

. Real gas version of the new Lagrangian method 

. Godunov scheme and TVD scheme 
. Numerical examples for supersonic real gas flow 

. Concluding Remark 


New Lagrangian formulation and the conservation 
form for 2-D steady supersonic flow (perfect gas) 

Independent variables of Euleriau formulation: x , y 

Independent variables of the new Lagrangian formulation: r , £ 
(Hui and Van Roessel, 1985) 

t Lagrangian time 

£ stream function 


Lagrangian conservation forin:(Loli and Hui, 

<9E OF n 
+ “ 


J.C.P., V.89, 207, 1990) 

( 1 ) 
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where 


/ 


E = 


V 


X \ 

0 ^ 

H 


0 

Ku + pV 

F — 

—pv 

Kv-pU 

5 * — 

pu 

U 


—u 

v y 


\ -v / 


IL — 


u = 


dx 

dr 

dx 


v. = 


V 


dy_ 

dr 

dy 

dt 


K = p(uV - vU ) 

H — (u 2 + v 2 )/2 + h(p,p ) 


h(p,p) = 


IP 


(7 ~ 1 )P 


Advantages of the new Lagrangian method 

( 1) Truly multi-dimensional computation (no time splitting), a com- 
putational cell is literally a fluid particle, physics is closeldly followed. 

( 2) Crisp slip line (contact ) resolution^ 1 point) 

( 3) No grid generation is needed, grids are automatically generated 
following streamlines. 

( 4) Particularly appropriate for hypersonic computation 
( 5) Ready for parallel computation. 


Real gas version of 
the new Lagrangian method 

. Accurate real gas Riemann solver 
. Decoding of E to find u,v,p,p 
Tannehill’s EOS is used ( NASA CR 2470, 1974) 
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initialize 



u.v.p.p, 
U arnlV 



Godunov/ TV D scheme 


Rieinann solver 


compute 
flux vecor F 


time lines 

, V, / 


streamlines 




Physical plane 


Computational plan 


compute new 
vector E 


> 


upgrade 
E by TVD 

Cell j 



limiter (Sweby) 

zzu 



TVD scheme 

boundary 


decode E and 
output u.v.p. p, 

U and V 


output x and y 
coordinates 


Computational Domain and Mesh 


Pressure, density and temperature contours for the real gas chan 
nel problem; (a) isobars, (b) isopycnics 





























typical time line 



Pressure, density and temperature contours for a sophisticated 
real gas outlet problem, showing interactions among the waves and 
shock reflection on the body surface; (a) schematic sketch, (b) isobars, 
(c) jsopvnics a ml (d) isotherms. 


Concluding Remarks 


Accurate yet efficient real gas effect computation: 

As accurate as Tannehill’s EOS warrants 
Genuine steady approach saves CPU 

Bears all the advantages of the Lagrangian method mentioned above 
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A TIME-ACCURATE IMPLICIT METHOD FOR CHEMICAL 
NON-EQUILIBRIUM FLOWS AT ALL SPEEDS 

Jian-Shun Shuen 
Sverdrup Technology, Inc. 

Lewis Research Center Group 
NASA Lewis Research Center 


Summary 

A new time-accurate coupled solution procedure for solving the chemical non- 
equilibrium Navier-Stokes equations over a wide range of Mach numbers is de- 
scribed. The scheme is shown to be very efficient and robust for flows with velocities 
range from M < 10“ 10 to supersonic speeds. 

Description 

Chemically reacting flows in aeropropulsion systems are often not amenable 
to numerical algorithms developed for high speed compressible flows. Examples 
include the rocket motor flow which involves a wide range of Mach numbers (from 
near zero velocity at the closed end to supersonic at the nozzle exit) and flow in gas 
turbine combustor where velocity is in the low subsonic range, yet the variation in 
density is very large so as to preclude an incompressible approach. Unfortunately, 
most available solution methods are usually applicable to either compressible flow 
with moderate to high Mach numbers or incompressible flows, but not both. This 
difficulty points out the need for a single solution technique that can be employed 
in flows with wide Mach number range and large property variations. 

There are two underlying reasons for the difficulty of low Mach number compu- 
tations in compressible flows. When the Mach number becomes low, the eigenvalues 
of the system differ widely so that the equations becomes stiff, resulting in signifi- 
cant slowdown in the convergence rate. Another reason is the singular behavior of 
the pressure term in the momentum equations as Mach number approaches zero, 
causing large roundoff error and smearing the pressure variation field, and, conse- 
quently, preventing the convergence. 

In the present numerical algorithm, methods to overcome these difficulties are 
proposed and successfully tested. The resulting code is then extended to consider 
the effects of non-equilibrium chemistry and real-gas properties. Included in the 
abstract are results from a linear stability analyses for the present algorithm (refered 
to as the pressure based scheme) and a conventional compressible flow algorithm 
(density-based algorithms) , and convergence histories for the convergent-divergent 
nozzle flows. The transonic nozzle flow has an inlet Mach number of 0.03 and exit 
Mach number 3.0. The subsonic cases are obtained by adjusting the nozzle back 
pressure so that flow is not choked at the throat. These results demonstrate the 
efficiency of the present method for flows over a wide Mach number range. 
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MOTIVATION 


• Reacting flows in aeropropulsion devices are often not amenable to 
numerical algorithms developed for high speed compressible flows, e.g., 

- rocket motor - wide range of Mach numbers, from near zero velocity 
at closed end to supersonic flow at nozzle exit. 

- gas turbine combustor - low subsonic velocities, but large variation in 
density precludes incompressible approach. 

• Most popular codes for gas turbine combustor and other low-speed re- 
acting flows are TEACH-type codes, which were based on technologies 
developed more than 20 years ago. 

• Compressible flow CFD technologies can be carried over to low-speed 
flow regime. 


OBJECTIVES 

1. Development of an efficient and robust Navier-Stokes code for chemi- 
cally reacting flows at all speeds. 

2. Implementation of liquid-fuel spray combustion and turbulent combus- 
tion closure models. 

3. Demonstration of the capability of the code for aeropropulsion appli- 
cations, including flows in gas turbine combustor, after burner, liquid 
propellant rocket combustion chamber, etc. 
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STABILITY ANALYSIS. CFL~tOQ 


STABILITY ANALYSIS. CFL*100 



Magnitude of the Maximum Eigenvalue of the Amplification 
Matrix G, Pressure- Baaed Scheme, CFL=100. 




Magnitude of the Maximum Eigenvalue of the Amplification 
Matrix G, Density- Baaed Scheme, CFL=100. 



Convergence History of the Pressure-Based S cheme for the Convergence History of the Density-Based Scheme for 

C-D Nozzle Flow, CFL=10O. C-D Nozzle Flow, CFL=100. 
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Difficulties with Compressible Flow Algorithms at Low Mach 
Numbers 


• Disparities among system’s eigenvalues (stiffness), u, u + c, u — c, re- 
sulting in significant slowdown in convergence rate. 

• Singular behavior of pressure gradient term in momentum equations 
as Mach number approaches zero, 


P*u * 2 + 


P* 


As Mach number is decreased, pressure variation (Ap* oc M 2 ) becomes 
of similar magnitude as roundoff error of the large pressure gradient 
term 


METHOD OF APPROACH 

Pressure Singularity Problem 

• Pressure decomposed into two parts: 


P^Po + Pg 


p g replaces p in momentum equations and retains p g as one of the 
unknowns. 

• Employs conservative form of the governing equations, but uses prim- 
itive variables 

(Pg,U,vJl,Yi) 

as unknowns. Conservation property preserved for flows at moderate 
to high Mach numbers. Pressure field accurately resolved for low Mach 
number flows. 

Eigenvalue Stiffness Problem 

• Pressure rescaled so that all eigenvalues have the same order of mag- 
nitude. Physical acoustic waves removed and replaced by a set of 
pseudo-acoustic waves which travel at speed comparable to fluid con- 
vective velocity. 
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Discretized Equations 

{r - ArD + + Ar(|* - ^)) AQ = -Ar(R)» 


where 


3QP-4CT + Q"- 1 3(E-E V ) 

” 2A t + dx 


-H 
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3E 

3Q’ 


_ 3H 
D = — 

dQ 


A - dE " 

-A* V ~ - ! 

dQ 

dQ 


T = 


3Q 


• Use dual-time stepping technique for time-accurate solutions. 

• Fully implicit, fully coupled solution method. Use LU factorization for 
multi-dimensional flows. 


ONE DIMENSIONAL GOVERNING EQUATIONS 
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Maximum Eigenvalue of 


CONVERGENCE HISTORY FOR THE C-D NOZZLE FLOW, CFL = 1()<). 


Non-Equilibrium Air Dissociation/R.ecombination Cliemistry 
5 Species, 11 Reaction Steps 

CD No'^le. Pressure Bosed Algorithm CD No^rie. Density Bosed Algorithm 




STABILITY ANALYSIS 

Linear stability analysis of 1-D inviscid equations using Euler implicit time 
inarching technique. 


STABILITY ANALYSIS. CFL=100 



Magnitude of the Maximum Eigenvalue of the Amplification 
Matrix G, Pressure- Based Scheme, CFL^IOO. 


stability analysis. cfl=ioo 



Magnitude of the Maximum Eigenvalue of the Amplification 
Matrix G, Density-Based Scheme, CFL=100. 
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CONVERGENT-DIVERGENT NOZZLE FLOW, Mi=0.01 


Non-Equilibrium Air Dissociation/Recombination Chemistry 
5 Species, 11 Reaction Steps 





COMBUSTION OF n-PENTANE FUEL DROPLET 


Non-Equilibrium n-Pentane - Air Combustion Chemistry 
5 Species, 1 Global Reaction Step 



a-a 



Motor Froclion 


DROPLET COMBUSTION, M ~ 10“ 4 - 10" 6 


Non-Equilibrium n-Pentane - Air Combustion Chemistry 
5 Species, 1 Global Reaction Step 



DROPLET COMBUSTION, M ~ 10~ 4 - 10 -6 
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CONCLUSIONS 


• An efficient and robust numerical algorithm has been developed for 
chemically reacting flows at all speeds. 


FUTURE PLAN 

1. Development and validation of the 2D/3D all-speed code for reacting 
flows. 

2. Implementation of zonal capability. 

3. Implementation of turbulence, turbulence/reaction closure, spray 
combustion, and thermal radiation models. 

4. Reacting flow simulation for gas turbine combustor, rocket motor 
combustion chamber, etc. 
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Study of Shock-Induced Combustion Using an Implicit TVD Scheme 


Shayne Yungster 

Institute for Computational Mechanics in Propulsion 
NASA Lewis Research Center 

The supersonic combustion flowfields associated with various hypersonic propul- 
sion systems, such as the the ram accelerator (a ramjet-in-tube concept), the oblique 
detonation wave engine, and the scramjet, are being investigated using a new CFD 
code. The code solves the fully coupled Reynolds-averaged Navier-Stokes equations 
and species continuity equations in an efficient manner. It employs an iterative 
method that is based on the lower-upper symmetric successive overrelaxation (LU- 
SSOR) implicit factorization scheme, and a second order symmetric total variation 
diminishing (TVD) differencing scheme. To accelerate the convergence of the basic 
iterative procedure, this code can be combined with vector extrapolation methods, 
such as the Minimal Polynomial (MPE) and the Reduced Rank (RRE) Extrapola- 
tion. The extrapolation procedure solves a linear least squares problem and produces 
a sequence of approximations that, in general, has better convergence properties than 
the sequence obtained from the iterative scheme alone. Two different formulations of 
the LU-SSOR factorization scheme are currently implemented. In one formulation, 
the implicit operator includes the full Jacobian matrix of the chemical source term, 
leading to a preconditioner matrix of size n t x n,, where n, is the number of species, 
which has to be inverted at every grid point. If the number of species considered 
is large, inverting this preconditioner can be very expensive. Therefore, a second 
formulation has been introduced in which the Jacobian matrix is replaced by a di- 
agonal matrix that is designed to approximate the time scaling effects obtained by 
using the full Jacobian. In this formulation, no matrix inversions are required. Fig- 
ure 1 shows the density residual history obtained with the two formulations (with and 
without extrapolation) for the case of a supersonic flow past a compression corner. 
The chemical nonequilibrium processes are simulated by means of a 9-species, 18-step 
finite-rate H 2 ~a. 1 T combustion model. When extrapolation is used, it is started after 
No iterations, and is implemented in the so called “cycling” mode, using a sequence of 
K m ax vectors obtained from the iterative scheme. The overhead in CPU time due to 
the use of extrapolation is very small (less than 1%) in the present case. The results 
indicate that savings of up tp 40% in the overall computational work required to 
reach convergence can be realized by using RRE in combination with the basic iter- 
ative scheme (using both formulations). Similar results are obtained with MPE. The 
results also indicate that for the present chemistry model, the diagonal formulation 
is less efficient, requiring approximately 45% more CPU time than the Full Jacobian 
formulation to reach convergence. The code is currently being applied to study shock 
wave/boundary layer interactions in premixed combustible gases, and to investigate 
the ram accelerator concept. Results obtained for a ram accelerator configuration 
(Fig. 2) indicate a new combustion mechanism in which a shock wave induces com- 
bustion in the boundary layer, which then propagates outwards and downstream. 
The combustion process creates a high pressure region over the back of the projectile 
resulting in a net positive thrust force. 
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Objective 


• To investigate shock- wave/boundary layer interactions 
involving premixed combustible gases. 


• To compute (based on the Navicr-Stokes equations) the 
supersonic combustion flowficld in a ram accelerator, a 
ramjet-in-tube concept. 


W 


(b) 




Figure X: Supersonic reacting flow past a compression corner; (a) pressure con- 
tours; .(b) Convergence History of L 2 density residual (grid, 80 x 50). 


Pressure 



cur 

Accelerator Barret 




Temp ere tore 


Figure 2: Nondimensional pressure (top) and temperature (bottom) contours 
for a ram accelerator, a ramjet-in-tube device for accelerating projectiles to ultrahigh 
velocities. Mixture: stoichiometric H 2 ~ air, M — 6.7, P « = 1 atm, T*, = 300° K, 
T w = 600° K. The vertical direction is magnified by a factor of 2. 
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Equations and Transport properties 

• 2D/axisymmctric Reynolds- Averaged Navier-Stokos equations 

(general coordinates). 

• Species specific heat, conductivity and viscocity obtained from fourth 

order polynomials of temperature. 

• Conductivity and viscosity of mixture - Wilke’s mixing rule. 

• Binary mass difiusivities - Chapman-Enskog theory. 

Combustion and Turbulence Models 

• 9-species, 18-st.ep chemistry model for hydrogen-oxygen combustion. 

• Bakhvin-Lomax Turbulence Model 

Numerical Method 

• Finite difference method. 

• LU-SSOR implicit factorization scheme. 

• Second-order symmetric TVD differencing scheme. 

• Use' of extrapolation techniques for convergence acceleration 

(i.e. Reduced Rank and Minimal Polynomial Extrapolation). 


Reduced Rank Extrapolation Method (RRE) 

• Given the vector sequence Q°, Q 1 , ..., Q i+1 , with k = Kmax , com- 
pute the differences 

AQ J = Q J+1 - Q j , j — 0, 1, k. (1) 

• Next, determine the scalars 70 , 71 ,. .-, 7 i- by solving the constrained 
least-squares problem 

minimize || E 7jAQ j '|| (2) 

j=<> 

subject to Ej_ 0 7 j = 1 . 

• Finally, set 

Siu- = £ 7 ,q j o) 

Although the original definition of RRE is different , the definition given 
above is equivalent to it and results in an implementation that is more 
stable numerically. 
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Numerical Method 


• LU-SSOR implicit factorization scheme: 

LT- l USQ =s -At[R/fS] 

L = I + $At[D^A + + D~B + - A~ — B~ — C\ 

U = I + /3At[Df A" + D+B- + A + + B + \ 

T = 1 + 0At[A + + B + — A~ — fl“] 

• Second-order symmetric TVD differencing scheme: 

[RHS\ = (F Mjt - F Htk + G ] k+i - + W jJc 

+ + R * h *^ - *h*h> 

Where 

*!f+J = -Q' + }) 

* W ” \ (z* + «*)/2«, if|a|<«i 

= *iitf i+4 i + ivji + 

.1 < 6 < .4 


Reacting flow past a compression corner 




96 





Convergence History of L 2 Density Residual 



Schematic of "Oblique Detonation” Ram Accelerator 

(L = 15 cm; d p = 2.5 cm; d t = 3.8 cm; 0 = 14°) 
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NONDIMENSIONAL PRESSURE (TOP) AND TEMPERATURE (BOTTOM) CONTOURS 
Stoichiometric Hydrogen-Air; M=6.7; P=1 atm; T=300 K. 


PRESSURE 
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TEMPERATURE 


NONDIMENSIONAL TEMPERATURE (TOP) AND MACH NUMBER (BOTTOM) CONTOURS 
Stoichiometric Hydrogen-Air; M=7.2; P=1 atm; T =300 K. 
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NONDIMENSIONAL TEMPERATURE (TOP) AND MACH NUMBER (BOTTOM) CONTOURS 
Stoichiometric Hydrogen-Air; M=8.0; P=1 ntm; T=300 K. 


TEMPERATURE 
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Conclusions 


A new CFD code that efficiently solves the Navier-Stokes equations 
with finite-rate chemistry was presented. 

The application of vector extrapolation methods to viscous, chemically 
reacting flows was demonstrated. 

Results indicate that significant savings in computational work can be 
realized by using vector extrapolation methods in combination with 
the basic iterative scheme. 

Ram Accelerator Concept: results indicate that viscous effects are 
of primary importance. The combustion processes in the boundary 
layer strongly affect the entire flowfield. 
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UPWIND SCHEMES AND BIFURCATING SOLUTIONS 
IN REAL GAS COMPUTATIONS 


Ambady Suresh 
Sverdrup Technology, Inc. 

Lewis Research Center Group 
NASA Lewis Research Center 
and 

Meng-Sing Liou 

Internal Fluid Mechanics Division 
NASA Lewis Research Center 

The area of high speed flow is seeing a renewal of interest due to advanced 
propulsion concepts such as the NASP, Space Shuttle, and future civil transport 
concepts. Upwind schemes to solve such flows have become increasingly popular in 
the last decade due mainly to their excellent shock capturing properties. Within 
this class of upwind schemes the Osher scheme has a few distinct advantages, such 
as a continuously differentiable numerical flux (making it especially attractive for 
implicit schemes), entropy satisfying solutions and the exact resolution of stationary 
shocks and contacts. 

High speed flow is generally accompanied by changes in the chemical and ther- 
modynamic composition of the gas which need to be taken into account. In the first 
part of this paper we present the extension of the Osher scheme to equilibrium and 
non-equilibrium gases. For the equilibrium case, the extension proceeds by choos- 
ing the pressure and entropy as the independent thermodynamic variables so that 
the Riemann invariants can be reduced to quadratures. Efficiency is achieved by 
noting that the quadrature need be only as accurate as the discretization error. In 
the non-equilibrium case, the determination of the intermediate points using Rie- 
mann invariants involves iterative solution of a transcendental equation. As this is 
computationally inefficient, an approximate method which involves no iteration is 
presented. For simplicity, the source terms are treated explicitly. In both cases, the 
formal extension is unique unlike the Roe scheme where there is a one parameter 
family of solutions with no clear choice among them. 

Computations based on the above scheme are presented to demonstrate the 
feasibility, accuracy and efficiency of the proposed scheme. One of the test problems 
is a Chapman- Jouguet detonation problem for which numerical solutions have been 
known to bifurcate into spurious weak detonation solutions on coarse grids. Our 
results indicate that the numerical solution obtained depends both on the upwinding 
scheme used and the limiter employed to obtain second order accuracy. For example, 
the Osher scheme gives the correct CJ solution when the super-bee limiter is used 
but gives the spurious solution when the Van Leer limiter is used. With the Roe 
scheme the spurious solution is obtained for all limiters. 



OVERVIEW 

• Extension of the Osher Scheme for reacting flows. 

• 2 species system 

• Intermediate and Sonic points 

• Numerical results 

• Bifurcating Z — N — D Detonations. 

• Theory 

• Numerical Results 


THEORY [Colella, Majda, Royturd (86)] 


Simplified Model: 


u t + (r w 2 + Ahz\ = 0 

\2 / if 


u~p/p 


Z n — k(j)(ll)z 


For large kAx and/or Ah non-physical weak detonation waves 
exist which travel at mesh speed 


NUMERICAL RESULTS 


• Riemann Problems (ideal gas) 

• Reacting nozzle flow 
•Z — N — D detonations 
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Soni c P o in t s 


N 



Using invariants is expensive due to exponents. 

r - M 

' (Ai-AJr) 


= %(1 - t.) + C/«r s 


Freeze eigenvectors at (Ul+Ur)/ 2 and evaluate wave strengths 
as: 


An,4 = ^[Ap^poAu] 

Aro = f-Ap + a 2 Ap] 
a 1 

Ar 3 = A z . 


Intermediate states obtained by transalation from left and 
right. 
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HRE 


N e 2 O e z P 



For multispecies Uy and Up can be obtained only iteratively: 

ap a N + bp# = c 

However, Uy, Up need only be as accurate as the overall dis- 
cretization error. 


1) Source terms do not change wave pattern. 



2) Along x = 0, solution approaches homogeneous 
Riemann problem near origin. Li and Yu [1985] 


- STRATEGY 


Use homogeneous Riemann problem [Gottlieb, etc.] for up- 
winding and add in source terms later. 
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X 


Interface flux F i+ i not well defined. 
However . . . 
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BIFURCAT ING Z-N-D DETONATIONS 


Problem 



r(z) — constant. 

Solved on moving grid. 
k - Step function of T. 


NUMERICAL RESULT S 

• Bifurcation depends on limiter used and Riemann Solver. 


R. S. \ Limiter 

I s * Order 

L.L. Limiter 

Superbee 

Osher-N 

B 

B 

NB 

Osher-R 

B 

B 

NB 

Roe 

B 

B 

B 
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PRESSURE OENSITY 

-M M Ltf 0.00 Q.2S 0.50 0.7S 1.00 


IDEAL GAS SHOCK TUBE PROBLEM 
OSHER SCHEME NATURAL ORDER 



Fig. 2: Ideal gas shock tube problem. Frozen solution. 
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OENSITY 


C - J DETONATION PROBLEM 
OS HER SCHEME NATURAL ORDER 







PRESSURE DENSITY 


C - J DETONATION PROBLEM 
OSHER SCHEME NATURAL ORDER 





Fig. 4a: Chapman-Jouguet detonation problem on a coarse grid. 
Reaction zone = 1/10 cell. Superbee Limiter. 
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PRESSURE DENSITY 


C - J DETONATION PROBLEM 
OSHER SCHEME NATURAL ORDER 



Fig. 4b: Chapman- Jouguet detonation problem on a coarse grid. 
Reaction zone = l/io cell, Vanleer Limiter. 
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DIVERGENT NOZZLE 

OSHER SCHEME NATURAL ORDER 



Fig. 5: Divergent Nozzle Problem. Reaction zone = 20 cells. 
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CONCLUSIONS 


• The Osher scheme can be extended efficiently to reacting 
flows. 


• Bifurcations in reacting flows require futlier study. 
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HIGH ORDER PARALLEL NUMERICAL SCHEMES 
FOR SOLVING INCOMPRESSIBLE FLOWS 


AviLin 

University of Pennsylvania 

Edward J. Milner, May-Fun Liou, and Richard A. Belch 
Internal Fluid Mechanics Division 
NASA Lewis Research Center 

The use of parallel computers for numerically solving flow 
fields has gained a lot of importance in recent years. In most 
cases, the parallel machine executes modified standard serial 
CFD codes, taking advantage of parallelism by using special 
constructs, directives, system commands and calls. In other 
instances, numerical schemes and computational algorithms have 
been changed or redesigned to increase performance of the 
overall algorithm. This paper presents a new high order 
numerical scheme for CFD specifically designed for parallel 
computational environments. 

A distributed MIMD system gives the flexibility of treating 
different elements of the governing equations with totally 
different numerical schemes in different regions of the flow 
field. This heterogeneous parallel numerical approach is 
sometimes technically complicated in terms of programming and 
administering the various processes, but it is quite efficient 
and contains enough free parameters for optimizing execution 
speedup . 

The parallel decomposition of the governing operator to be 
solved is the primary parallel split. The decomposition of the 
physical domain of the flow into subdomains is the secondary, or 
optional, split. At the coarsest parallel level, each of the 
processors is assigned to solve a suboperator of the original 
PDE operator over a given subdomain. An iterative numerical 
procedure is employed. All of these splittings are designed to 
ensure a high rate of convergence. 

The primary parallel split was studied using a hypercube-like 
architecture having clusters of shared-memory processors at each 
node. The approach is demonstrated using as examples some 
simple steady-state incompressible flows. Future studies should 
investigate the secondary split because, depending on the 
numerical scheme that each of the processors applies and the 
nature of the flow in the specific subdomain, it may be possible 
for a processor to seek a better, or higher order, scheme for 
its particular subcase. 
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OBJECTIVE: 


TO REDUCE THE TIME REQUIRED TO SOLVE 
CFD PROBLEMS BY USING PARALLEL 
PROCESSING TECHNIQUES 
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DISCUSSION FRAMEWORK 


• PARALLEL PROCESSING TECHNIQUE FOR 

SOLVING A BLOCK TRIDIAGONAL SYSTEM 

• SOME APPLICATIONS & RESULTS 

• PARALLEL PROCESSING TECHNIQUE FOR 

INVERTING A MATRIX 
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A Parallel, Scalable, 2D Navier-Stokes Solver 


• Parallel algorithm design provides efficient implementation on distributed 
or shared memory machines 

- 74% efficiency achieved on Hyper cluster test bed 

• Scalable and portable algorithm 

- basic building blocks ( e.g. matrix solvers ) 

- total solver 

• Second-order accurate 



Stream function, driven cavity, SO by 50 grid 
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ASSUME THAT 



B = K * B 


m,l m,l 


1.1 



B 


m,2 




B 


m,3 


K * B 

m,3 


3,3 
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THEN 

1) (A u + A! 2 *K 2fl + A, 3 * K 3 ,)* B| tl = I 

2) (A 2 , * K 1<2 + A 22 + A 2i3 *K 3 , 2 )* B 22 — I 

3 ) (A 31 * Kj 3 + A 3|2 * K 23 + A 3 3 ) * 63,3 — I 

4) A| t Aj 2 + A 1>3 * K 3i2 — 0 

5) A 2 1 * Kj 3 + A 2 2 A 2 3 — 0 

6) A 31 + A 3 2 "ft JC 2 j + A, 3 *©~ 0 

7) Au * ^1,3 Aj 2 * K 2 .3 + Aj 3 = 0 

8) A 2 1 + A 2 2 * K 2i + A 2 3 * K 3 1 = 0 

9) A 3 1* Kj 2 + A 3 2 + A 3 3 * K 3>2 — 0 

COMPUTATIONAL TECHNOLOGIES BRANCH 

= 1MASA =====^^ LEWIS RESEARCH CENTER: 

THUS, 

1) 1,, = (A w + A,, 2 *g] + A, 3 *(g))-‘ 

2) i 2 . 2 = (a 2i1 *(g)+ a 2 , 2 + a 2 , 3 * PQ )- 1 

3 ) I3.3 = (A 3 , * H + A 3 , 2 * (g)+ A 3 , 3 )-‘ 

4 ) @=-A 1 i ,,*(A,, 2 + A u *[jQ) 

3) — A 2 2* (A 2 ] * K 13 + A 23 ) 

6 ) (g)= -A‘.3 * (A 3 , + A 3 2 * H ) 

7) K i3 = (A, j — A, 2 * A 22 * A 2 *(A 12 * A 22 *A 2/ 3 — Aj 3) 

8) ^2,1 = (A?„2 “ A 2 3 * A3.3 * Aj 2 ) '* (A 2 3 * Aj 1 3 * A 3 , — A 2 ,) 

9 ) |K 3 , 2 h (A3.3 - Aj,, * A*, * A li3 y i * (A 31 * Aj,* A, 2 - A 3 , 2 ) 
COMPUTATIONAL TECHNOLOGIES BRANCH^ 
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SUMMARY 


THE PARALLEL SOLVER: 
+ + + ): 




1) IS GENERAL 

2) CONVERGES VERY QUICKLY 

3) IS READILY EXPANDABLE 

a) CAN BE USED WITH ONLY ONE PROCESSOR 

b) WILL USE AS FEW OR AS MANY PROCESSORS 
AS ARE AVAILABLE 


1) PERFORMS SOME REDUNDANT CALCULATION 
BECAUSE OF OVERLAP OF THE DATA 
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BLOCK TRIDIAGONAL SYSTEM 
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BLOCK TRIDIAGONAL SYSTEM 
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BLOCK TRIDIAGONAL SYSTEM 
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HYPERCLUSTER TEST BED ARCHITECTURE 
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FLOW NEAR A ROTATING DISK 


EQUATIONS: 


&u i *p . J dhl . 0 ( u \jL dlu \ 

-T + "to = -7& : + M'S ? + a : W + ^ r ) 

, ur , Bv /««. 8 / v\ , S*.\ 

+ T + U, 8r = y \8^ + IFl l Tj + &r} 


<ho . 

«ST + * 


ITERATION 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 


1 dp , \ dho , 1 dxo , 3ho 

z ~Jdz^~ V \dr^r dr + 

i + + = 0. 

• ‘ r * Sz 


I 


RESIDUAL ERROR 

MAXIMUM ERROR - 0.15019E+02 
MAXIMUM ERROR - 0.52285E+01 
MAXIMUM ERROR - 0.23837E+01 
MAXIMUM ERROR - 0.10327E+01 
MAXIMUM ERROR - 0.39308E+00 
MAXIMUM ERROR - 0.94527E-01 
MAXIMUM ERROR - 0.47686E-02 
MAXIMUM ERROR - 0.14872E-04 
MAXIMUM ERROR - 0.78692E-09 
MAXIMUM ERROR = 0.60396E-13 



Flow In Ihe neighborhood of s disk 
rotating in a fluid at rest 

vaocre cw ewfnt K w-radw. v-drewm w a. » 1 
A lay«f of Mt) It can<«d by tM M4ng to M acS 
or vneewt tore**. TM ctrttotofal tort*# to »** Wn It 
glvt rtH to ttcondtry HowwWcft It dhcM nrStaty 
outward. 
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FLOW NEAR A ROTATING DISK 


• DISTRIBUTED MEMORY SIMULATION 

• SERIAL MATRIX INVERTER 




Flow in the neighborhood of a dbk 
rotating in a fluid at reft 

Vtocfy ca wpfw n l fc u-mm, v-Otuw« i >«n»* » , *-wul 

AlayworfcMtecvrtotytwOMiomtittgtMacftxi 

•T tore**. Th* etnMtftf tetM to Itt Mn lay*r 

gMtrtMtoMcoMttytowwhtehkiarMMrwMy 
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i v# uj/ y 

BOUNDARY LAYER ALONG A PLATE 

EQUATIONS: 

du , du Bhi 

Uco 

u ei + v ai = v W‘ 

u a y\ 

®3i _i_ — o 

0* + By U * 

y =0: u == v = 0 ; y —oo: 

^ ^°° The boundary layer along a flat 

plate at zero incidence 

ITERATION 

RESIDUAL ERROR 

1 

MAXIMUM ERROR - 0.14258E+01 

2 

MAXIMUM ERROR - 0.26890E+00 

3 

MAXIMUM ERROR - 0.28632E-01 

4 

MAXIMUM ERROR - 0.24873E-03 

5 

MAXIMUM ERROR = 0.16154E-07 

6 

MAXIMUM ERROR > 0.11084E-11 
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MATRIX INVERSE 
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A 
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B 

B 
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3.1 

3.2 

3.3 


3.1 

3.2 

3.3 
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BOUNDARY LAYER ALONG A PLATE 


• DISTRIBUTED MEMORY SIMULATION 

• SERIAL MATRIX INVERTER 



NUMBER OF PROCESSORS 



The boundary layer along a flat 
plate at zero incidence 
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A Parallel, Scalable, 2D Navier-Stokes Solver 



Stream function, driven cavity, 50 by 50 grid 

• Parallel algorithm design provides efficient implementation on distributed 
or shared memory machines 

- 74% efficiency achieved on Hyper cluster test bed 

• Scalable and portable algorithm 

- basic building blocks ( e.g. matrix solvers ) 

- total solver 

• Second-order accurate 

— COMPUTATIONAL TECHNOLOGIES BRANCH " 
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RELATIONSHIPS 


A * B + A * B + A * B = I 



1.1 
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u 
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34 
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An Efficient and Robust Algorithm For 
Time Dependent Viscous Incompressible Navier-Stokes Equations 


John W. Goodrich 

Computational Fluid Dynamics Branch 
NASA Lewis Research Center 


A recently developed finite difference algorithm is presented for unsteady in- 
compressible Navier-Stokes calculations. The algorithm is extremely robust with 
respect to Reynolds number, and has been used to directly compute incompressible 
flows with smoothly resolved streamfunction, kinetic energy and vorticity contours 
for Reynolds numbers as high as Re ~ 100,000 without requiring any subscale 
modelling. The accompanying figure shows contour plots for the streamfunction 
and vorticity at nondimensional times t — 83 and t = 84 from a transient calcu- 
lation at Re = 25,000, with a 256 X 256 grid and At = The algorithm is 

second order accurate in both time and space, with a Crank-Nicolson differencing 
for the diffusion terms, with a lagged second order Adams-Basforth differencing for 
the convection terms, and with central differencing for all space derivatives. There 
is no constraint on the time step size from diffusion time scales, but the convection 
time scales impose a stability constraint of — < 1, and typically a CFL number of 
0.75 or 0.80 is used. This algorithm is based on the fourth order Partial Differential 
Equation fof incompressible fluid flow which uses the streamfunction as the only 
dependent variable, 


d A ip 

dt 


dip dip dip dtp 1 2 

+ -^-A-—- - —A—- - —A 2 ip = 0, 
dy dx dx dy Re 


for x in 0, and t > 0. 


Notice that the vorticity does not enter into this formulation. Although the al- 
gorithm currently is only for incompressible flows in two space dimensions, the 
algorithm can be and is being extended to handle temperature dependant density, 
subsonic compressible flows, and flows in three space dimensions. The algorithm 
produces discretely divergence free velocity fields on a nonstaggered grid. The 
streamfunction discretization is related to a primitive variable discretization on a 
staggered grid, so that primitive variable velocity boundary conditions can be im- 
plemented for the discrete streamfunction. The algorithm is extremely efficient with 
respect to use of both CPU time and physical memory. A typical time dependant 
calculation with a 128 x 128 spatial grid and At - — • requires approximately 2 
MegaBytes of memory, and approximately 2.3 CPU seconds per time step on an 
IBM RS/6000 model 530 workstation. Codes implementing the algorithm have been 
run for problems of various scales on systems ranging from PC’s and workstations 
to CRAY supercomputers. The algorithm has been implemented on an INTEL 
TOUCHSTONE parallel processor. Solutions will be shown for cavity and channel 
flows at various Reynolds numbers. 
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TIME DEPENDENT INCOMPRESSIBLE NAVIER-STOKES EQUATIONS 

— + V • (uu) -Au = -Vp + F, for x in Q, and t > 0, 

dt ' Re 

V • u = 0, for x in 0, and t > 0. 


STREAMFUNCTION EQUATION FOR UNSTEADY INCOMPRESSIBLE FLOW 


with 


d A ip X . , , dip .dtp dip dip , . n , . 

— - — A 7 ip+ -^-A-~ - for x in 0, and t > 0; 

dt Re dx dy dy dx 

u(x,t) = and v{x,t) = for x in fl, and t > 0. 

oy ox 


* Initial and boundary conditions must be supplied. 


=* For n C R 3 . 

=> Vorticity and pressure DO NOT enter into the streamfunction formulation. 
==> A single equation for a single scalar unknown. 

=> The velocity solution is always divergence free, and so incompressible. 


THE STREAMFUNCTION ALGORITHM FOR UNSTEADY INCOMPRESSIBLE FLOW 




-HP) + £«<*•)-•£ 


At 

+ T 


S x [S y (z n ) L&(z n )J-6 v [6 x (z'')U(z’') 
6. ( (5„(£ n-1 )La(z n_1 )'] - fu*" -1 ) La(£- 1 ) 


with 


2Ay 


«y+i and <y = 


2Ax 


K 


+ 1.J 




* z n « t/>(*,t n ) on the discrete grid. 

* The discretization uses central spatial differencing throughout. 

* All variables are defined at each grid point. 

* Banded LU decomposition and multigrid solvers have been used. 


=> The variable xp is as smooth “as possible”. 

=> The implicit equation is elliptic, for all Reynolds numbers. 

=> The local domain of dependance is the large symmetric 13 point discretization stencil. 
=> The discrete solution is exactly incompressible, <5 x (u ( "V) + S y (v"V) = 0. 

=> The stability limit is a CFL constraint, < 1- 
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A MULTIGRID SOLVER FOR THE LINEAR IMPLICIT EQUATIONS 


La(z" + 1 ) — Bi(£ n + 1 ) = Source Ttrm{i n , 1 ) 

Ail 1C 

* The Biharmonic operator is factored as two Laplacians. 

* Point Gauss-Seidel smoothing (or Red Black Gauss-Seidel). 

* Linear restriction and prolongation. 

* A V-cycle with 3 iterations per grid level while coarsening, none while refining. 

=> uj — A 0 is introduced ONLY for iterating with the biharmonic operator. 
==> Implemented on scalar, vector and parallel systems. 


TYPICAL PERFORMANCE DATA 

CAVITY PROBLEM (IBM RS/6000 model 530 workstation) 

* Square cavity transient from t = 0 to t = 1, at Re = 9600. 

* 10 to 15 iteration cycles reduce residuals to less than 5.0 x 10’ 13 . 

(1) 256 x 256 grid, A t = 7 grid levels, 6.6 MBytes storage, 10.2 scc/ts; 

(2) 192 x 192 grid, A t = 7 grid levels, 3.8 MBytes storage, 5.5 sec/ts; 

(3) 128 x 128 grid, A t = 6 grid levels, 1.7 MBytes storage, 2.5 sec/ts. 

=> CPU time increases “linearly” with the number of grid points: 

(1) 1.54 x 10" 4 sec/ts per grid point on the 256 X 256 grid; 

(2) 1.48 x 10“ 4 sec/ts per grid point on the 192 x 192 grid; 

(3) 1.50 x 10“ 4 sec/ts per grid point on the 128 X 128 grid. 

=> Required memory increases “linearly” with the number of grid points: 

(1) 104.78 Bytes per grid point on the 256 x 256 grid; 

(2) 106.97 Bytes per grid point on the 192 x 192 grid; 

(3) 107.12 Bytes per grid point on the 128 x 128 grid. 

=> The number of iteration cycles/ts is independent of resolution and Reynolds number. 
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STREAMFUNCTION MAX AND MIN COMPARATIVE DATA 
The square driven cavity at Re=5000 


The Streamfunction Minimum 


Source 

Grid 

'/’min 

*r» 


Vm trt 

Kim and Moin 

96 x 96 

-1.12 x 10" 1 




Goodrich 

128 x 128 

-1.15 x 10“ 1 

66 _ 
128 

132 

260 

60 _ 
128 

Goodrich 

256 x 256 

-1.18 x 10 _1 

132 

260 


137 

260 

Ghia, Ghia, and Shin 

256 x 256 

-1.19 x lO" 1 

131 

260 


137 

260 


The Streamfunction Maximum 





^Pm ax 

ax 


Urn ax 

Goodrich 

128 x 128 

3.44 x 10- 3 

102 _ 
128 

204 

260 

9 

128 

Goodrich 

256 x 256 

3.13 x lO' 3 

206 

266 


1 0 
256 

Ghia, Ghia, and Shin 

256 x 256 

3.08 x lO" 3 

207 

266 


1 9 
266 


TRANSIENT DRIVEN CAVITY 
Re = 25,000, 256 x 256 grid, t = 10.25 
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TYPICAL PERFORMANCE DATA 


CHANNEL PROBLEM (IBM RS/6000 model 530 workstation) 

* Backward facing half height step with a parabolic inflow profile on the step corner. 

* Transient from t - 0 to t = 100, Re = 800, u mai = 1.5 for inflow profile, At = 

* Outflow Boundary Conditions are: 


d 7 rl> d 7 w 

— - = 0 and — — = 0 O 
dx 7 dx 7 


0 and 


d -l± = 0 

dx‘ 


* 10 to 35 iteration cycles reduce residuals to less than 5.0 X 10“ 13 

* 30 iterations per cycle on the coarse grid. 

* Coded for 1024 x 32 grid with 4 grid levels using 3.5 MBytes storage. 

(1) 800 x 32 grid, 25 diameter domain, 4.7 s/ts; 

(2) 480 x 32 grid, 15 diameter domain, 3.1 s/ts; 

(3) 224 x 32 grid, 7 diameter domain, 1.7 s/ts. 


==> Relative CPU time decreases with the length of the domain: 

(1) 1.79 x 10" 4 sec/ts per grid point on the 800 x 32 grid; 

(2) 1.95 x 10" 4 sec/ts per grid point on the 480 x 32 grid; 

(3) 2.23 x 10" 4 sec/ts per grid point on the 224 x 32 grid. 


138 



STREAM FUNCTION CONTOURS 
Re=25000. 256*256 grid. 1 = 03.00 



STREAM FUNCTION CONTOURS 
Re=25000. 256*256 grid. t=04.OO 
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15 


V0RT1CITY CONTOURS 
Rc = 25000. '256*256 grid, t = 83.00 



VORTICITY CONTOURS 
Re = 25000. 256*256 grid. 1=64.00 
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o 5 ' o * 0 

I* 0 >7.0 
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STREAM FUNCTION CONTOURS 
Re=800, 64*448 grid, 0<=x<=7, t=30, ul=1.5 
Backward facing step with parabolic inflow, code g4a32vc3.f 
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BACKWARD FACING STEP 
Re = 2000, 32 x 480 grid, 0 < x < 15 
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STRE AMFUN CTION CONTOURS 
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SUMMARY 


* Vorticity and pressure DO NOT enter into the streamfunction formulation. 

* A single equation for a single scalar unknown in R 3 . 

* The variable ip is as smooth tt as possible”. 

=► Second order accuracy in both time and space. 

=> The discrete solution is exactly incompressible, <5 x (u?\) + S y (v”\) = 0. 

=}► The stability limit is a CFL constraint, < 1- 

==> Implemented on scalar, vector and parallel systems. 

=> CPU time increases “linearly” with the number of grid points. 

==> Required memory increases “linearly” with the number of grid points. 

=*► The number of iteration cycles/ts is independent of resolution and Reynolds number. 
=> Extremely robust with respect to Reynolds number. 


GOALS 

* Develop and adapt numerical methods. 

* Investigate dynamic phenomena (instable, transitional, unsteady time asymptotic). 
STEPS 

* Incompressible flow (robust, fast, validated). 

* Artifical outflow boundary conditions (Tom Hagstrom, Univ. of N.M.). 

* Parallel processor implementation on INTEL iPSC (Rodger Dyson, NASA Lewis). 

* Temerature dependant density (Boussinesq approximation). 

* Chemistry models (Marty Rabinowitz, NASA Lewis). 

* Multilevel adaptive methods (Steve McCormick, Univ. Co. Denver). 

* Compressible subsonic flow. 

* Higher order differencing scheme. 

* Three dimensional version. 

APPLICATIONS 

* Development of models for simplified chemical kinetics. 

* Development of models for chemistry /turbulence interaction. 

* Prediction of chemical emmisions (High Speed Research Program). 

* Understanding instabilities (flashback, flameout and accoustic/combustion interaction). 
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ADVANCES IN ENGINEERING TURBULENCE MODELING 


T.-H. Shih 

Institute for Computational Mechanics in Propulsion 
NASA Lewis Research Center 


Some new developments in two-equation models and second order closure mod- 
els will be presented. Two-equation models (e.g., k-e model) have been widely used 
in CFD for engineering problems. In most low- Reynolds number two-equation mod- 
els, some wall-distance damping functions are used, especially in the eddy viscosity, 
to acount for the effect of wall on turbulence. However, this often causes the con- 
fusions and difficulties in computing flows with complex geometry and also needs 
an ad hoc treatment neax the separation and reattachment points. In this paper, 
modified two-equation models axe proposed to remove abovementioned shortcom- 
ings. The calculations using various two-equation models axe compaxed with direct 
numerical simulations of channel flows and flat boundary layers. 

Development of second order closure model will be also discussed with empha- 
sis on the modeling of pressure related correlation terms and dissipation rates in 
the second moment equations. All the existing models poorly predict the normal 
stresses neax the wall and fail to predict the 3 dimensional effect of mean flow on the 
turbulence (e.g. , decrease in the shear stress caused by the cross flow in the bound- 
ary layer). The newly developed second order neax- wall turbulence model to be 
described in this paper is capable of capturing the near- wall behavior of turbulence 
as well as the effect of three dimension mean flow on the turbulence. 
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• Two-Equation Near- Wall Turbulence Models 

• Second-Order Closure Models 

• Modeling of Three-Dimensional Turbulent Flows 


k-e model 

The eddy viscosity ut is assumed in two-equation models as follows: 

_ , k 2 

v T = 

vt = Cpf^kr 


where r = k/e. 

The general k-e (or k-r) model equations are of the following forms: 


k,t H“ Ujkj — 

€ ,t == 

Tt + UjTj = 


V T 

Ok 

ut 

ut 


cr r2 


+ Kj 

+ A 

- + A T . 


j ,3 


+ II + ut (Uij + Uj ti ) — e + D 


+ (Uij + U,j) Uij - C,f 2 j + E 

+ r ( — + kjTj - - ( — + T.iT , 


J >7 


k \ a. 


T 1 


r \a r i 


+ (i - c tl )-MUi,j + + (c<2/ 2 - 1) 
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The damping function / M : 

/ M = /(^) = /(i? t ) — Jones and Launder model 

fp = /(-^) = f{y + ) — for all other k-e models 

The parameter R t and y + are used for counting the wall effect. 

The problems with the form of / M (y + ) are following: 

Near the separation region, — * 0, =>■ vt = 0; 

In the flows with complex geometry, the y is not well defined. 


y + is an unacceptable parameter. 



Turbulent kinetic energy ( Re 0 = 1410) 
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In the present study, we suggest another parameter: 

R= *H\v\ 

6 V 

and the following modifications in Shih’s model: 

ffi = 1 - exp |c 3 [l - exp(C 6 -R 1/4 )j } 



e = e [l - exp(-i? t 1/2 )j 

where R t = k 2 /ue, Cq = .004, C 3 = .0004, C 6 = 1.2 
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Turbulent flow past a hill 

(Re~1.4 . 10 6 ) 
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Another alternative: 


R u 


U 4 


ve 


and 


fti = 1 “ exp(-a 1J R^ /4 - a 2 Rl/ 2 - a 3 R u ) 

\CoVr k . 

e = e [l - exp(- J R t 1/2 )] 

where i? f = fc 2 /i/e, ax = 5 • 10 -3 , a 2 = TO -5 , a 3 = 8 • 10~ 7 , C 0 = .01. 


Second order modeling of near-wall turbulence 

The exact equation for the Reynolds stress tensor is: 

-j^(uiUj) = Pij + T{j + d\j + H-ij — £ij 

where ( ) stands for an ensemble average, D/Dt = d/dt + Ukd/dxk ■ 

Pij — {'U’iUk)Pj,k (’Uj'Ufc) Ui,k 

T{j = {UilljUk) ,k 

D if = v{UiUj) M 
n ij = --(uiPj + ujp ti ) 

r 

£ij — 2l 
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Using the near-wall asymptotic behavior of turbulence^ 10 ! as model 
constraints, we formed a set of modeled transport equations for the 
Reynolds-stress tensor and the dissipation rate of turbulent kinetic energy. 
The proposed near-wall model for II, • j — Eij is: 


II,y E{j 



[2(uiUj) +4((uiUk)njn k + (ujU k )nink)+2(ukUi)n k nininj] 


where n, is a unit vector normal to the surface, and f w — exp(— (i? t /Ci) 2 ), 

R t = Ci = 1.358_R 0 ^ 44 , R er = u t 6/u. u t is the friction velocity, 6 is 
the thickness of the boundary layer or the half width of the channel. 


Away from the wall, the velocity pressure-gradient correlation II; j is 
split into the rapid part Ely ^ and the slow part Ily^: 

n« = n!j> + nS? 

The proposed model for 11^ — is: 

- £ ij = -e(/% + - f w ) 

where 

P = 2 + + 80 ' 1 M 1 + 62.4(— // + 2.3 ///)]} exp(--^) 

F= 1 + 271 1 1 + 9 II 

III = hijbjkbm 

bjj = (uiUj)/((j 2 ) - 5ij/3 
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The rapid part of velocity pressure-gradient, 11^ is modeled as fol- 
lows(Shih and Lumleyl 11,12 !): 


nS)> =(g + 2a 5 K 9 2 )(i7 jJ - + Vj't) - |(1 - a s )(P,j - \piij) 

+ (I + T<iKAf - \Ptii) + ^Pa ~ Da) + | bijP 


+ 


where, 


3 3 

2 


Sit! 2 ) ( U j U k)Ui t q)(UkUq) ( UiUp)(UjUq)(Up t q + Uq t p}] 


Pij — (ni‘Ufc)J/j > fc {UjUk)Ui t k 
Dij = -(uiU k )U k j - {ujU k )Uk,i 

„ In 

P - 2 Pu 

as = ~i5 (1 + C2F ' /2) 

C 2 = 0.8[1 - exp(-(H t /40) 2 )] 


Dissipation rate equation 

The modeled dissipation rate equation derived in this work is: 

66 

e, t + Ui€ t i = (ue yi - - V»o^ 2 y 

« / 2 \ 

_ (q 2 ) ( UiU i^i,j ~ ^2 ~ {ukUi)(Uiji — Ui,ij)Ui t jk 

where 

= — + 0.98[l - 0.33 ln(l - 55//)]exp(-2.83^" 1/2 ) 

5 

iJji ~ 2.1 

= -.15(1 - F) 

y , 

4<<? 2 > 
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U, V , w 


U - MEAN VELOCITY 


K - KINETIC ENERGY 





Three-dimensional boundary layer 
The effect of three-dimensional flow on: 

• Reduction in Reynolds stress. 

• Lag development between the stress and strain rate. 
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- 0.4 -02 


Turbulence modeling: Launder and Shima model 


Turbulent shear stress Turbulent intensity 



Turbulence modeling: Lai and So model 


<uv> PROFILE 


u, v, w — RMS OF FLUCTUATING VELOCITY 



y+ 



y+ 
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Turbulence modeling: Present model 



Performance of existing models in modeling of 3-D flows: 

• Two-equation (e.g. k-e ) models — incapabale 

• Reynolds-stress models — possible 

Launder and Shima model — no 
Lai and So model — no 
Present model — yes 
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On Recontamination and Directional-Bias Problems 
in Monte Carlo Simulation of PDF Turbulence Models 
Andrew T. Hsu 
Sverdrup Technology, Inc. 

Lewis Reseach Center Group 
NASA Lewis Research Center 


Turbulent combustion can not be simulated adequately by conventional moment closure turbulence models. 
The difficulty lies in the fact that the reaction rate is in general an exponential function of the temperature, and 
the higher order correlations in the conventional moment closure models of the chemical source term can not be 
neglected, making the applications of such models impractical. The probability density function (pdf) method 
offers an attractive alternative: in a pdf model, the chemical source terms are closed and do not require additional 
models. 

The partial differential equation for the probability density function, f\ can be written as 

v 

pd t P + pv*d a P + 

»=i 
AT AT 

= -a a (p< v"h Pi> P) - pYm d li^(< ««!** > h 

jm 1 

where the terms represent the time derivative, mean convection, chemical reaction, turbulent convection, and 
molecular mixing, respectively. The fact that the pdf equation has a very large dimensionality renders finite 
difference schemes extremely demanding on computer memory and CPU time and thus impractical, if not entirely 
impossible. A logical alternative is the Monte Carlo scheme, wherein the number of computer operations increases 
only linearly with the increase of number of independent variables, as compared to the exponential increase in a 
conventional finite difference scheme. 

A grid dependent Monte Carlo scheme following that of J.Y. Chen and W. Kollmann has been studied in the 
present work. In dealing with the convection and diffusion of the pdf, the pdf equation is discretized on a given 
grid, e.g., 

P i x+irj = QjP rj>l + 0jP tj + 7 jP r,;-l 

where 

Oj + P) + 7; = 1 

However, if this is the only restriction satisfied by the numerical algorithm, the mass fractions may not be con- 
served due to re-contamination, and directional-bias also appears. These phenomena are illustrated in Figure 1: 
Consider a mixing layer; use white balls to represent contaminants in the upper stream and black balls to represent 
contaminants in the lower stream. As the two streams move toward right, the location of the white balls and black 
balls are interchanged randomly to simulate convection and diffusion. From Figure 1, it is clear that directional- 
bias caused by recontamination caused the center of the mixing layer to drift downward. ( Directional- bias can 
be partially corrected by changing sweeping directions. ) One also notices that after the first marching step, the 
conservation law is violated, reflected in the Figure as missing white or black balls. 

It is found that in order to conserve the mass fractions absolutely, one needs to add further restriction to the 
scheme, namely 

Qj + 7; = Qj-l + 7>+t 

A new algorithm was devised that satisfies this restriction in the case of pure diffusion or uniform flow problems. 
Using the same example, it is shown that absolute conservation can be achieved. This result is shown in Figure 2. 
One can see that the diffusion process is symmetric, and the problem of directional-bias is eliminated. 

Although for non-uniform flows absolute conservation seems impossible, the present scheme has reduced the 
error considerably. 
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CLOSURE PROBLEM: 


- Turbulent** Modeling 

- Analogy of shear stress: diffusion model. 
??? 

p id j j > • • • » ^ n j 2 , p') 

But in general: 

pwi ± pw i Y u ... ) Y, n T,p) 
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The Simulation of Convection and Diffusion 

wap = -oa(> < tew > h 

Discretize: 

r-\ T 

When' the coellicionts must satisfy: 

<*j + Pj + 7/ = 1 


MONTE CARLO METHOD 

(J ) Use an ensemble oi stochastic particles to represent the PDh 

(2) The particles move in such a way that the PDF evolution equa- 
tion is satisliod. 

(3) Three steps: eouveefion/dilhisiou; molecular mixing; chemical 
reaction. 
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PDF EVOLUTION EQUATION 


+ ftvji.P + (> 'ii 

#==1 

< r n\'h > P) - P 


D,;-' { Wj{ j/’i , 


;V .Y 
\ * V /)“ 

/=] 7 = 1 , ' ,r > 


-ph)P) 

(< <-ij\A- > P) 


(J) Moan Convection 

(2) Chemical Reaction 

(3) Turbulent Convection 
(•1) Molecular Mixing 


Dimension: N+5 


finite difference not feasible. 


PROBABILITY DENSITY FUNCTION MODEL 



ASSUMED PDF 


167 



MOMENT CLOSURE MODEL 


Example: 


Wf — —BT n exp ( — 7^7 ) p 2 YpYo 


T = T + T' 
Yf = Yf + Yf 
Y f = Y f + Yjr 


• w F can be expended into series and terms such as Y f Yq , YpT', 
YqT and T ' 2 have to be modeled. 

Restrictions: T' <C T, T„ ~ T 




1 


] + i 


j + 2 
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CONCLUDING REMARKS 


(1) For absolute conservation, the number of particles received in a 

cell must equal to the number that are sent out from the same 
cell: 

a j + 7 j = a j-i + 7j+i 

(2) The particles that are replaced must be the particles that were 

sent out. 

(3) Store at least 3 arrays of random indices and move them succes- 

sively. 

(4) Eliminated the problem of recontamination and directional-bias 

in the cases of uniform flow or pure turbulent diffusion. Reduced 
error. 

(5) In case of non-uniform flow, convection is needed, and thus (l) 

can not be satisfied. 

(6) The capability and accuracy of the present scheme are demon- 

strated through the example of a turbulent non-premixed flame. 
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Implementation of a k-e Turbulence Model to RPLUS3D Code 


Tawit Chitsomboon 

Institute for Computational Mechanics in Propulsion 
NASA Lewis Research Center 

The RPPLUS3D code has been, developed at the NASA Lewis Research 
Center to support the national aerospace plane project. The code has the 
capability to solve three dimensional flowfields with finite rate combustion 
of hydrogen and air. The combustion processes of the hydrogen-air system 
are simulated by an 18-reaction path, 8-species chemical kinetic mechanism. 
The code uses a Lower-Upper (LU) decomposition numerical algorithm as 
its basis, making it a very efficient and robust code. Except for the Jacobian 
matrix for the implicit chemistry source terms, there is no inversion of a 
matrix even though a fully implicit numerical algorithm is used. 

A k-epsilon turbulence model has recently been incorporated into the RPLUS3D 
code. Initial validations have been conducted for a flow over a flat plate. 
Results of the validation studies will be shown. Some difficulties in imple- 
menting the k-epsilon equations to the RPLUS3D code will also be discussed. 



MOTIVATION 


• MOST APPLICATIONS ESPECIALLY FLOWS IN SCRAM- 
JET COMBUSTORS ARE TURBULENT 

• COMPLEX FLOWS IN SCRAMJET COMBUSTORS (e.g., NOR- 
MAL H 2 INJECTION INTO SUPERSONIC AIR) WILL RE- 
QUIRE AT THE MINIMUM A k-e MODEL 


OBJECTIVE 


• IMPLEMENT A HIGH REYNOLDS NUMBER k-e TURBU- 
LENCE MODEL 

• A LOW REYNOLDS NUMBER MODEL IS NOT PRACTI- 
CAL FOR A 3D CHEMICALLY REACTING CFD CODE 


SOME FEATURES OF RPLUS3D 


• 3D CFD CODE DEVELOPED AT NASA LEWIS RESEARCH 
CENTER 

• FULLY IMPLICIT FINITE VOLUME L-U SCHEME 

. CENTRAL-DIFFERENCE SCHEME WITH JAMESON-TYPE 
A.V. 

• FINITE-RATE HYDRO GEN- AIR COMBUSTION (8-SPECIES, 
18-STEP) 

• EFFICIENT: NO INVERSION OF MATRICES EXCEPT FOR 
COMBUSTING FLOWS 

• DEVELOPED ESPECIALLY FOR NASP COMBUSTORS 

• COULD BE USED FOR GENERAL PROBLEMS BY SET- 
TING VARIOUS SWITCHES 



K-EPSILON TURBULENCE MODEL 

• HIGH REYNOLDS NUMBER FORM 

• ECONOMICAL FOR 3D COMBUSTING FLOWS 

• LAW OF THE WALL 

• SOLVE 2 PDE’S FOR k AND e 

• USE L-U ALGORITHM ( THE SAME AS THE FLOW SOLVER) 

• DECOUPLE FROM THE FLOW SOLVER 

• MODULAR PROGRAM : ONLY ONE SUBROUTINE CALL 
FROM THE FLOW SOLVER 


HIGH REYNOLDS NUMBER k-e MODEL 


• STANDARD k-e MODEL 

dpk dpkui 




a ' a* =W>-^v (, + S )v* 

d£e + dfXUi _ c pi _ C * + V • (/» + — )V £ 
at axi k k a e 


WHERE 


P = — UiU 


dui 

dxj 


AND, 


-UiUj = u t (uij + Ujj) - ( 2 / 3 )8ij(k + v t V • V) 
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BOUNDARY CONDITIONS FOR k-e MODEL 

• ASSUMPTIONS 

- LOG LAW PROFILE 

u + = ( l/K)ln{Ey+ ) 

- LOCAL EQUILIBRIUM 

dissipation — production 


• PSEUDO CODE 

- INTEGRATE FLOW EQUATIONS 

- INTEGRATE k-e EQUATIONS 

- FIND u T FROM LOG-LAW PROFILE BY NEWTON’S IT- 
ERATION 

- FIND k and e AT THE FIRST CELL CENTER 

k = (u T ) 2 /v^r : e = {u T ) 3 /ny 

-FIND v t = CJrh 


DIFFICULTIES ENCOUNTERED 

• RUN AT LOWER CFL NO. THAN LAMINAR CASES 

- DUE TO EXPLICIT VISCOUS TERMS IN THE FLOW 
SOLVER 

• IMPLICIT FORMULATION FOR BOTH CONVECTION AND 
DIFFUSION TERMS 

• PARTIALLY IMPLICIT SOURCE TERMS 

• POSITIVITY OF k-e IMPOSED AUTOMATICALLY 

e = e n k n+l /k n 

• REQUIRES ARTIFICIAL DAMPING 

• NO DAMPING AT TWO POINTS ADJACENT TO THE WALL 
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LAMINAR FLOW OVER A FLAT PLATE 


• MACH NO. = 0.5 ; 31 X 31 GRID 
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Figure . Turbulent boundary layer flow over a flat plate. 
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FUTURE PLANS 


. MORE VALIDATIONS 

- DOUBLE NORMAL JETS IN CROSS FLOW (UVA) 

- SLANT JET IN CROSS FLOW (NASA LEWIS) 

. COMPRESSIBILITY EFFECTS 
. PROBABILITY DENSITY FUNCTION (PDF) 
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TURBULENCE AND DETERMINISTIC CHAOS 


Robert G. Deissler 
LeRC Research Academy 
NASA Lewis Research Center 


ABSTRACT 

Several Turbulent and nonturbulent solutions of the Navier-Stokes equations 
are obtained. The unaveraged equations are used numerically in conjunction 
with tools and concepts from nonlinear dynamics, including time series, phase 
portraits, Poincare sections, largest Liapunov exponents, power spectra, and 
strange attractors. 

Initially neighboring solutions for a low-Reynolds-number fully developed 
turbulence are compared. The turbulence, which is fully resolved, is sus- 
tained by a nonrandom time- independent external force. The solutions, on the 
average, separate exponentially with time, having a positive Liapunov expo- 
nent. Thus, the turbulence is characterized as chaotic. 

In a search for solutions which contrast with the turbulent ones, the Reynolds 
number (or strength of the forcing) is reduced. Several qualitatively 
different flows are noted. These are, respectively, fully chaotic, complex 
periodic (see the accompanying figure), weakly chaotic, simple periodic, and 
fixed-point. Of these, we classify only the fully chaotic flows as turbulent. 
Those flows have both a positive Liapunov exponent and Poincar6 sections 
without pattern. By contrast, the weakly chaotic flows, although having 
positive Liapunov exponents, have some pattern in their Poincare sections. 

The fixed-point and periodic flows are nonturbulent, since turbulence, as 
generally understood, is both time-dependent and aperiodic. 

Besides the sustained (forced) flows, a flow which decays as it becomes 
turbulent is examined. As in the sustained case, the flow is extremely 
sensitive to small changes in initial conditions. The sensitivity increases 
with improved spatial resolution. For the finest grid (128 3 points) the 
spatial resolution appears to be quite good. 

As a final note, the variation of the velocity-derivative skewness of a 
Navier-Stokes flow as the Reynolds number goes toward zero is calculated 
numerically. The value of the skewness, which has been somewhat controver- 
sial, is shown to become small at low Reynolds numbers, in agreement with 
intuitive arguments that nonlinear terms should be negligible. 
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TURBULENCE PROBABLY STILL MAIN UNSOLVED PROBLEM IN FLUID 
DYNAMICS. 

BUT CONSIDERABLE ADVANCES BEING MADE: 

o NUMERICAL SOLUTIONS OF INSTANTANEOUS 
NAVIER-STOKES EQUATIONS 

o UTILIZATION OF CONCEPTS AND TOOLS FROM NONLINEAR 
DYNAMICS, INCLUDING DETERMINISTIC CHAOS. 

"DETERMINISTIC" REFERS TO THE EVOLUTION EQUATIONS 
(NO RANDOM COEFFICIENTS.) 

THESE TWO GO HAND IN HAND. IN FACT, USE OF NONLINEAR DYNAMICS IN STUDY 
OF TURBULENCE (OR IN THE STUDY OF ANYTHING) DID NOT GET FAR UNTIL ADVENT 
OF HIGH-SPEED COMPUTERS, BECAUSE THE NONLINEAR EVOLUTION EQUATIONS 
(NAVIER-STOKES EQUATIONS IN CASE OF TURBULENCE) COULD NOT PREVIOUSLY BE 
EFFECTIVELY USED. 


HERE WE WANT TO TAKE A LOOK AT TURBULENCE BY USING CFD, TOGETHER WITH 
CONCEPTS AND TOOLS FROM NONLINEAR DYNAMICS. 


OBTAINED SEVERAL SOLUTIONS OF INSTANTANEOUS NAVIER-STOKES 
EQUATIONS AT LOW REYNOLDS NUMBERS-BOTH TURBULENT AND NONTURBULENT 


FLOW KEPT FROM DYING OUT BY TIME-INDEPENDENT FORCING TERM 
ALTHOUGH FORCING WAS STEADY, RESULTING FLOW WAS OFTEN 
HIGHLY UNSTEADY AND RANDOM IN APPEARANCE 


THIS IS A REMARKABLE FEATURE OF MANY TURBULENT FLOWS - E.G., A 
TURBULENT (TIME DEPENDENT) FULLY DEVELOPED PIPE FLOW MAY BE 
PRODUCED BY A STEADY PRESSURE GRADIENT 
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SOLVED NUMERICALLY THE UNAVERAGED THREE-DIMENSIONAL 
NAVIER-STOKES EQUATIONS WITH A TIME- INDEPENDENT 
FORCING TERM 

SUj- aiuiUfel .! ap ^ 6 2 Uj 

dt dx k p dXj dx^ 1 

3 

WHERE Fj = £ " a i V cos 0 ' Q • X 

n=l 


WHERE 


Vk(2 f l,l) 2 a r k(l,2,l) 3 a r k(l,l,2) 1/2 

Re a = ( u 1 ) x 0 /v 

2qj - (1,-l.D 3 qj=(U,-l> 


AND 


i a 2 p . y^ic 1 

p dXjdXj d X|X k 


PLOT OF FORCING TERM Fj ON PLANE 
THROUGH GRID CENTER 


MAGNITUDE OF FORCING VECTOR ON PLANE 
THROUGH GRID CENTER 
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USED A CUBICAL COMPUTATIONAL GRID (0 TO 2n ON A SIDE) 
WITH FOURTH-ORDER SPATIAL DIFFERENCING. 32 3 GRID POINTS 


USED A PREDICTOR-CORRECTOR TIME DIFFERENCING (SECOND-ORDER LEAPFROG 
PREDICTOR AND THIRD-ORDER ADAMS-MOULTON CORRECTOR) 

BOUNDARY CONDITIONS: PERIODIC 


OBTAINED ASYMPTOTIC (LONG-TERM) SOLUTIONS 


TO INSURE EXCITATION OF AS MANY MODES AS POSSIBLE 
AT A GIVEN REYNOLDS NUMBER, INITIAL CONDITIONS ON 
MOST OF ‘FLOWS WERE SPATIALLY CHAOTIC 


MAGNITUDE OF CHAOTIC INITIAL VELOCITY VECTOR 
ON PLANE THROUGH GRID CENTER 


|U| 



2k 
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CALCULATED SPATIAL VARIATION OF VELOCITY 
ON PLANES THROUGH GRID CENTER 



PLOT OF PROJECTION OF VELOCITY-VECTOR FIELD 
ON x r x 2 PLANE THROUGH GRID CENTER 
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Time series 


ASYMPTOTIC 
REYNOLDS NUMBER 





SIMPLE PERIODIC FLOW 
•2 r— Re a = 6.24 

JVTTVXVTTVTTV' 

— | I*" Period 

*.6 1 l 1 1 L I I 1 I 

90 95 100 105 110 115 120 125 130 


Re- = 6.72 




128 138 148 158 168 170 

Time 


Phase portrait 


U2(rc,rc,n) 
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R# a « 6.60 



u^onm^ijrfie.aan/ie) 
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U 2 {9re/8. 21)1/16, 23n/16) 


COMPLEX PERIODIC FLOW 



= 6.93 


1.0 1.2 


WEAKLY CHAOTIC FLOW 


FULLY CHAOTIC FLOW 







SUMMARY 


FIXED POINT 

LIMIT CYCLE, SIMPLE PERIODIC (PERIOD 1) 

LIMIT CYCLE, COMPLEX PERIODIC (PERIOD 2) 

WEAKLY CHAOTIC FLOW (LIAPUNOV EXPONENT POSITIVE, POINCARE' SECTION HAS 
SOME PATTERN) 

LIMIT CYCLE, COMPLEX PERIODIC (PERIOD 4) (PERIODIC WINDOW) 

FULLY CHAOTIC FLOW (LIAPUNOV EXPONENT POSITIVE, POINCARE' SECTION HAS NO 
APPARENT PATTERN, TURBULENT) 
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REYNOLDS NUMBER 







INCOMPRESSIBLE TAYLOR VORTICES 


If the angular velocity of the inner 
cylinder is larger than a critical 
value, Couette flow becomes unstable 
so that a secondary flow starts to 
appear with nonvanishing radial 
and axial velocity components; 
this new flow, called Taylor-Couette 
flow, is in the form of opposite 
toroidal vortices arranged next to 
each other in the axial direction. 


Axisymmetric laminar 
Taylor vortices 



COMPRESSIBLE TAYLOR VORTICES 


In the presence of heat transfer and angular motions, the influence of compress- 
ibility is nontrivially coupled and its effect on stability has yet to be determined. 


The addition of compressibility to the stability analysis greatly increases its com- 
plexity, since even in the linearized analysis there will now exist fluctuations 
not only in velocity and pressure but also in density, temperature, viscosity and 
thermal-conductivity. 
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GOVERNING EQUATIONS 


• Mean Flow 

momentum: 


energy: 


OP to? 

Or r 


sKfc-7)h?(S : ?)- 


• Boundary Conditions 


at r = l: w = 1, T - 1, P = /» = 1 

andat r-_: w = —,T=-=i> 


GOVERNING EQUATIONS 


• Small Perturbation 

<f>U, V,(> t) = <s> + <j>(ii)cxp[i{a£ +PC- wt)] 


where 

£, ??, £: computational coordinates in axial, radial 
and circumerential directions, respectively 
<£: steady unperturbed flow, 

— [t/^, X* 0 , /i^, , k 0 ] 

<j>\ complex disturbances 
a: disturbance wave number in £ direction 
0: disturbance wave number in C direction 
u: complex frequency 
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# Boundary Conditions for Disturbances 

- at r = 1 : 

1 i 

and at r = : 

Rl 

Note that no boundary conditions are required for pressure fluctuations when the 
equations are solved on a pressure staggered mesh. 


u, 5, u>, f = 0 
u, h, w, f = 0 


METHOD OF SOLUTION 


• Numerical solutions for the mean flow 

® Chebyshev collocation spectral method for stability analysis 
Momentum and Energy Equations: 

Gauss-Lobatto points for u,v,vj,t 
Continuity Equation: 

Gauss points for p 


Momentum and Energy Equations: 

Agl + B gi ,I/gz, (<t> 4 I g L? ) + 4 I g l P) 

= w [d G i,L G £,(^ 4 lg L P) 4- 4 


Continuity Equation: 

B G L G lg^ 4 C G (lg^ 4 P) = u>B G (lg^ 4 P) 


where 


/u\ 

v 

w 5 

\ f / 


p - (p) 


and 

A, B, C, D, E: coefficient matrices 

L: spectral differentiation operator 
I: spectral interpolation matrix 
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Velocity 


VERIFICATION 


Cojjipnrisou of the Critical Reynolds Number for Temporal Stability of Taylor - 
Couette Flow (i2J/i7J=0.5) 


Sparrow el ul™ Present 


n;/n; 

0.0000 

0.0600 

0.1250 

0.1800 

0.2350 

- 0.10 

- 0.20 

-0.30 

-0.35 

-0.50 


Re c 

68.19 

73.02 
84.30 

107.21 

221.33 

66.18 

70.03 
80.93 
88.81 

114.24 


Rc c 

68.187 

73.422 

34.346 

100.923 

220.867 

66.172 

69.895 

80.418 

89.947 

115.025 


RESULT AND DISCUSSION 


MEAN VELOCITY PROFILES 


!?J//tj=0.5 and H;=0 


(a)T;/T;=0.5 


(b)T;/T;=i.o 


WS/TSml* 
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Tempcralure 


RESULT AND DISCUSSION 


MEAN TEMPERATURE PROFILES 


/?;//?;= 0.5 and q-=o 




(b)Tj - /T;«1.0 


(c)r,vr;=i.5 


U»lO - 01 

W»0 5 

tf-f.O 

W-JO 

if-5 0 












RESULT AND DISCUSSION 


EFFECT OF TEMPERATURE RATIO 


<1=3.1 02, n;=u 






RESULT AND DISCUSSION 


EFFECT OF SPIN MOTION 


i\[ =1.0 and ife=500 


MAXIMUM GROVmi RATE 
VS 

TEMPERATURE RATIO 
a=6.2, Af= 1.0, and 72c=5 00 
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RESULT AND DISCUSSION 


m — i .o, Rc=m), n;/nj=-o.5 

o=G.2, atul i?=1.0 


Radial- Velocity fluctuation 



CONCLUDING REMARKS 


• For the case of rotating inner cylinder, the higher Mach number flows are unstable 
to a larger range of wave numbers. 

• By increasing the temperature ratio, the flow between counter rotating cylinders 
is stabilized with the growth rate decaying in a linear fashion. However, the effect 
of increasing temperature ratio is to amplify the maximum growth rate in the flow 
with a positive rotating speed ratio. 

• Further examinations can be made on those effects of axial flow, gap width, finite 
axial length, and other possible influential factors. 
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Techniques for Animation of CFD Results 

Dr. Jay Horowitz 
Computer Science Division 
NASA Lewis Research Center 
and 

Jeffery C. Hanson 
Sverdrup Technology, Inc. 

Lewis Research Center Group 
NASA Lewis Research Center 

Intro 

Video animation is becoming increasing vital to the CFD researcher, not just 
for the presentation results, but for recording and comparing dynamic visualizations 
that are beyond the current capabilities of even the most powerful graphic worksta- 
tions. To meet these need Lewis Research Center has recently established a facility 
to provide users with easy access to advanced video animation capabilities. How- 
ever, producing animation that is both visually effective and scientifically accurate 
involves various technological and aesthetic considerations that must be understood 
both by the researcher and by those supporting the visualization process. Many of 
these considerations will be addressed in the presentation. 

Scan Conversion 

Conversion of high-resolution workstation images to low- resolution television 
images is performed by hardware that either converts only a portion of the screen, 
or invokes various pixel-averaging techniques. The former can result in excessive 
aliasing, whereas the latter can result in smearing and disappearance of points and 
lines. Both can have trouble displaying semi-transparent objects from workstation 
using bit-mask transparency algorithms. Scan conversion is a major problem in 
accurately portraying Computational Fluid Dynamics grid geometry. 

Color Conversion 

Translation of the workstation’s red-green-blue component color space to the 
composite NTSC color system can generate various artifacts. This can be a major 
problem when color is used to convey scientific information. Several techniques are 
used at LeRC overcome these problems including: color-table adjustment to avoid 
over-saturation, and the use of edge-outlined contour plots to emphasize subtle color 
differences. 


Spatial Ambiguities 

Since passive viewers are denied the interactivity of workstations, that allows 
researchers to accurately determine spatial relationships and resolve visual ambigu- 
ities, it is critical that recorded data be easily understood. Careful attention must 
be paid to relative object and ’camera’ movements and scene content. In addition, 
advanced display techniques such as solid, smooth-shaded stream-tubes instead of 
line-segment tracers can significantly improve flow visualization. 
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Video as a CFD Research Tool 


* Visualize Time-Dependent Results 

* Record Interactions with Data 

* Workstations too slow for complex data 

* Video cost effective & ubiquitous 


TV Monitor 





Successful CFD Video Visualizations: 

ledmokigX’ 

• Scan Conversion 

• Color Conversion 

• Recording media 

Aesthetics - 

• Resolving Spatial Ambiguties 

• Motion Control 


Problems: 


^ me mmvare&f 

steps necessary to go from workstation 
to video 

W- * '*:?■' ?: : ;V/ • •• 

Most visualization soft 
provide all necessary 


m 


• I 

* ’ Video-ifying ’ the visualization usually 


Scanline 


Red Signal 


, /s//s//s/sss/sss//ss//sssssss/s//s//sssA 

Green Signal 


' ?///////////> 



Blue Signal 


Workstation RGB Color Domain 








Conversion 




j i .uuj i mi— p . - 


r* Video 


Hi-res Workstation, 


Scanline 


.uniinancc 


Chrominance Signal 


Saturation 


Composite Signal 

NTSC Color Domain 

Susceptible to adjacent pixel color interference 
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Effect of Scan Conversion 
on Texture Patterns 



NTSC Color Artifacts 

* Zippers 

* Dot Crawl 

* Saturation 


Image at risk 
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Modify Colormaps to avoid over-saturation 


Luminance: 

Y = 0.29*r + 0.59*g + 0J1 *b 
Chrominance: 

I = 0.59*r - 0.27*g - 032*b 
Q=0.21 *r - 0.52*g + 0 .31 *b 
Saturation: 

S = (I 2 +Q 2 ) 112 


O.K. If -0.2S < (Y-S) , (Y+S) < 1 .0 
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advanced Graphics & Visualization laboratory’s 
digital video facility 


RGB Monitor 
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Techniques for CFD Animation 

* Use thick lines when possible 

* Outline contour edges if possible 

* Use nonsaturating colors 

* Use shaded tracers 

* Preview images on TV monitor 

* Avoid simultaneous viewpoint and 
data motion 



DISTRIBUTED VISUALIZATION FOR COMPUTATIONAL FLUID DYNAMICS 


Don J. Sosoka and Anthony A. Facca 
Computer Science Division 
NASA Lewis Research Center 


Distributed concurrent visualization and computation in Computational Fluid 
Dynamics is not a new concept. Specialized applications such as RIP 
(Realtime Interactive Particle-tracer) and vendor specific tools like DGL 
(Distributed Graphics Language) have been in use for some time. Little work 
however, has been done to put together a concise, easy to use set of routines 
which would allow the CFD code developer to take advantage of this 
technology. Recent advances in the speeds of the supercomputer and graphics 
workstations have made it possible to realize many of the applications for 
distributed visualization which have long been envisioned. Among these are 
the ability to "watch" a CFD code as it runs on a supercomputer, to perform 
interactive grid generation on a local workstation interacting with a solver 
running on the supercomputer, and ultimately, to "steer the CFD run 
interactively from the workstation. For any of these applications to be 
developed, a set of tools are required to provide this distributed capability to the 
scientist. Three components can be identified as being necessary to make this 
capability useful to the CFD researcher. First, a FORTRAN-callable set of 
network routines must be provided to "link-up" the distributed processes and 
provide continuous process-to-process communication. Data is moved from the 
computational process to the visualization process transparently. Secondly, a 
FORTRAN-callable set of visualization functions must also be provided to the 
CFD researcher. These functions must allow the researcher to easily view the 
data in a meaningful, understandable way. Finally, some empirical data should 
be available to provide guidance to the researcher on the type of applications 
which would benefit from this kind of distributed processing. Several problems 
must be addressed in providing each of these three components. In dealing 
with the network, capacity is always a concern. Binary data compatibility 
between machines with vastly different architectures is a major factor which 
must be addressed. Most importantly, ease of use from the CFD code 
developers point of view must be maintained in delivering a useful solution. 

This paper describes a current project underway at NASA Lewis Research 
Center to provide the CFD researcher with an easy method for incorporating 
distributed processing concepts into program development. Details on the 
FORTRAN callable interface to a set of network and the visualization functions 
are presented along with some results from initial CFD case studies that employ 
these techniques. 
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OBJECTIVE: 


PROVIDE CFD RESEARCHERS AN EASY METHOD TO 
INCORPORATE DISTRIBUTED VISUALIZATION TECHNIQUES 
INTO THEIR PROGRAM DEVELOPMENT. 


DISTRIBUTED VISUALIZATION: 

CONCURRENT GRAPHICAL RENDERING ON ONE SYSTEM 
OF RESULTS, AS THEY ARE COMPUTED ON ANOTHER. 


Distributing Computational Code 


EB BB 


Cray Y/MP 



SGI 4D/340 
! 
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APPLICATIONS FOR DISTRIBUTED VISUALIZATION IN CFD 

• Convergence ' Status Checking (PEEKER) 

• Interactive Batch Grid Generation 

• Immediate Visualization of Time-Dependent Calculations 

• Interactive "Steering" of Remote Computations 


Supercomputer Personal Workstation 




DEGREES OF DISTRIBUTED VISUALIZATION 

• All computation, scientific visualization, and graphics 
performed remotely. 

ex: Ultra Frame Buffer 

• Computations and scientific visualization done remotely, 
graphics performed by the workstation. 

ex: DGL for SGI workstations 

• Numeric Computations done remotely, all scientific 
visualization and graphics performed by the workstation. 

most difficult to implement 

• requires network interface capable of 
transferring transparent binary data 

highest payoff 

• provides greatest degree of parallelism 

• no static datasets required to visualize results 
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TOOLS REQUIRED FOR DISTRIBUTED VISUALIZATION 


• FORTRAN-Callable Set of Network Interface Routines 

• establish network connection between processes 

• provide continuous process-to-process communication 

• FORTRAN-Callable Set of Visualization Functions 

• incorporate state-of-the-art rendering technology 

• utilize non-proprietary graphics interface 

• encompass all current CFD post-processing capability 

Design Goal: Make Routines Easy to Use From a CFD Code 
Developers Point of View 


Open Close Routines: 

OPNCON -- Open a "connect" on a port 
OPNLIS - Open a "listen" on a port 


Send Routines: 

PUTRA1 -- Send array of real numbers 

PUTIA1 -- Send array of integer numbers 

PUTSTR -- Send character string 

PUTRA2 -- Send 2-dimensional array real numbers 

PUTIA2 -- Send 2-dimensional array of integers 


Receive Routines: 


GETRA1 

GETIA1 

GETSTR 

GETRA2 

GETIA2 


Receive array of real numbers 
Receive array of integer numbers 
Receive character string 
Receive 2-dimensional array real numbers 
Receive 2-dimensional array of integers 



c Sample SEND program. 


c establish network connection 
CALL OPNCON ( ••• ) 



CALL PUTRA2 ( • • • ) 

• 

• 

Network 

Interface 


• 






c Sample RECEIVE program. . 

• 

0 

0 

Network 

Interface 



c accept network connection 
CALL OPNLIS ( ••• ) 


c receive 2D array from server 
CALL GETRA2 ( • • • ) 


XDR Decoder 

T 



c *** SVP/NETWORK RECEIVE CODE 

PARAMETER ( ID=8 0 , JD=70 ) 

REAL xgrd(ID, JD) ,ygrd(ID, JD) 

REAL ws(ID,JD) 

c *** access the network 

call OPNLIS 

I c *** get data from the network 
j 100 CONTINUE 

call GETIAl(l,nx, istat) 
call GETIA1 (1 ,ny , istat ) 
call GETRA2 (nx , ny , 2 , ws , istat) 
call GETRA2 (nx,ny,2, xgrid, istat) 
call GETRA2 (nx , ny , 2 , ygrid, istat ) 

c *** introduce data to SVP 



} * 

call GDATA (nx,ny,nx, ny, xgrid,ygrid) i 

call FDATA (ws, ' ') 


call GENCON ( ' ') 

call CLRSVP 
GOTO 100 

STOP 

END 
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TECHNIQUES FOR GRID MANIPULATION AND ADAPTATION 


Yung K. Choo 

Internal Fluid Mechanics Division 
NASA Lewis Research Center 

Peter R. Eisemann 
Program Development Corporation 
and 

Ki D. Lee 

University of Illinois 


Two approaches have been taken to provide systematic grid manipulation for 
improved grid quality. 

One is the control point form (CPF) of algebraic grid generation. It provides 
explicit control of the physical grid shape and grid spacing through the move- 
ment of the control points. The control point array, called a control net, is a 
sparse grid-type framework in physical space. The CPF-based grid manipula- 
tion has many merits. The method is efficient and easy to implement, since its 
formulation is algebraic and concise. Grid quality can be easily enhanced by 
numerous local and global grid distribution strategies. It is also compatible with 
various complementary operations needed in a quality grid-generation proce- 
dure. It works well in the interactive computer graphics environment and hence 
can be a good candidate for integration with other emerging technologies. 

The other approach is grid adaptation using a numerical mapping between the 
physical space and a parametric space. Grid adaptation is achieved by 
modifying the mapping functions through the effects of grid control sources. The 
source strengths are extracted from the distribution of flow or geometric 
properties on the initial grid. Grids can be made adaptive to geometry, flow 
solution, or grid quality, depending on what property is used to define the 
source strengths. One advantage of the method is that the basic characteristics 
of the initial grid can be retained while adapting it to grid quality. The use of grid 
control sources allows for linear combinations of different controls based on the 
superposition principle. And the grid can be adapted to more than one property 
through a series of mappings. The source formulation promotes smooth 
variations in the grid, even with irregularly distributed sources. The adaptation 
process can be repeated in a cyclic manner if satisfactory results are not 
achieved after a single application. 
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OUTLINE 


1. EXPLICIT GRID MANIPULATION BY USER 

TECHNICAL APPROACH - CONTROL POINT FORM (CPF) 
OBJECTIVES / ADVANTAGES 
STATUS /RESULTS 
FUTURE DIRECTION 

2. GRID ADAPTATION TO SOLUTIONS, QUALITY 

TECHNICAL APPROACH - NUMERICAL MAPPING 
with GRID CONTROL SOURCE 

OBJECTIVES / ADVANTAGES 

STATUS / RESULTS 

FUTURE DIRECTION 


OBJECTIVES 

TO DEVELOP EASY, EXPLICIT GRID MANIPULATION CAPABILITY 


ADVANTAGES 

THE CPF-BASED GRID GENERATION: 

- IS EFFICIENT AND EASY TO IMPLEMENT 

- CAN EASILY ENHANCE GRID QUALITY BY VARIOUS LOCAL/GLOBAL 
GRID DISTRIBUTION STRATEGIES 

- CAN BE USED AS A COMPLEMENTARY TOOL AS WELL AS A STAND-ALONE TOOL 

- WORKS WELL IN THE INTERACTIVE COMPUTER GRAPHICS ENVIRONMENT 

- ALLOWS TO MANIPULATE GEOMETRY (FREE-FORM BOUNDARY) 

- MAKES INTERACTIVE CFD ATTRACTIVE 
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TECHNICAL APPROACH - CONTROL POINT FORM (CPF) 
CONSTRUCTION OF A CURVE & A CONTROL NET 



Q ( r ' t ) = T(r, t) + a,[l - Gj(r)] [P(l,t) - Fj(t)J 
+ a 2G«_i(r)^P(N - l,t) - F N (t ) j 
♦«3p -Hi(t)][P(r,l) -ERr)] 
+ «4 H H-i(t)[P(r,M - 1) - E H (r)] 


TURBO/I INTERACTIVE PROCESS 





GEOMETRY MODIFICATION USING FREE-FORM BOUNDARY 





TURBO/I AS A COMPLEMENTARY TOOL 
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FUTURE DIRECTION 


SURFACE GRID CONTROL USING CONTROL POINT FORM 
INTERACTIVE ADAPTIVE GRID GENERATION 

TECHNICAL APPROACH 
NUMERICAL MAPPING WITH GRID CONTROL SOURCE 




(s,t) parametlic domain 
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Grid control sources definition: 



where = weighting factors 

<{> = a grid quality parameter 



where (s’ , f ’) = modified parametric coordinates 

K s ljkl ,K l l j kl = influence coefficients 
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OBJECTIVES 


TO DEVELOP GRID ADAPTATION CAPABILITY TO SOLUTIONS, GRID QUALITY, 
AND GEOMETRY TO OBTATIN ACCURATE PREDICTION OF COMPLEX FLOWS 


ADVANTAGES 


THE APPROACH IS EASY TO IMPLEMENT AS ONE OF THE REDISTRIBUTION 
SCHEME. 

THE USE OF GRID CONTROL SOURCES ALLOWS TO ADAPT THE GRID BY USING 
LINEAR COMBINATIONS OF FLOW PROPERTIES. 

THE APPROACH CONTROLS BOTH GRID DENSITY AND QUALITY WITHOUT 
ANY ADVERSE INFLUENCE OF ONE TO THE OTHER. 

THE GRID ADAPTATION SCHEME CAN BE INTEGRATED WITH WELL-DEVELOPED 
FLOW SOLVERS AS SUBROUTINES. 


NEEDS 

FULL 3D ADAPTATION WITH REDUCED COST 


STATUS / RESULTS 

2D ADAPTATION TO SOLUTIONS, GRID QUALITY 
APPLICATION OF 2D ADAPTATION TECHNIQUE TO 3D PROBLEMS 
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CORNER FLOW 
GRID COMPARISON 



Adapted grid 
@ x = 0.407 


Adapted grid 
@x = 0.688 

CORNER FLOW 
SOLUTION COMPARISON 


Adapted grid 
@ x = 1.000 



Density contours 
on initial grid 
@ x=* 0.407 


Density contours 
on initial grid 
@ x= 0.688 


Density contours 
on initial grid 
@ x = 1.000 
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FUTURE DIRECTION 


EXTEND IT TO FULL 3D ADAPTATION 


EXPLORE AND IMPLEMENT COST REDUCTION SCHEMES 


EXPLORE POSSIBLE BENEFITS FROM PARALLEL COMPUTING 


APPLY IT TO CHALLENGING PROPULSION PROBLEMS 


INTERACTIVE GRID GENERATION 

Program TURBO 



CONSTRUCT A SIMPLE CONTROL NET 
GENERATE SURFACE GRID 


COMPUTE & EXAMINE INITIAL 
GRID 



SELECT A CONTROL POINT TO 
MODIFY GRID 




EXAMINE NEW GRID COMPUTE NEW GRID 

-ACCEPT OR REPEAT THE PROCESS 



TRANSLATE THE CONTROL POINT 
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